

ABSTRACT 

The purpose of the study was to demonstrate the capability of predicting 
two-dimensional, compressible and reacting flow in the combustion chamber and 
nozzle of the Space Shuttle Main Engine (SSME). CHAM's general purpose 
Computational Fluid Dyanmics (CFD) code, PHOENICS, has been used. A non- 
orthogonal body fitted coordinate system has been used to represent the nozzl 
geometry. The Navier-Stokes equations are solved for the entire nozzle with 
the k^e turbulence model. The wall boundary conditions have been calculated 
based on the wall functions which account for pressure gradients. 

Results of the demonstration test case reveal all expected features of the 
transonic nozzle flows. 

Of particular interest is the location of an internal shock, and regions of 
highest temperature gradients. Calculated performance (global) parameters 
such as thrust chamber flow rate, thrust and. specific impulse are also in 
reasonable agreement with available data. 

Recommendations are made for further improvements in the physical models, as 
well as for the use of numerical model at NASA MSFC. 
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Section 1 


INTRODUCTION 

1 . 1 Thrust Chamber Simulation Status 

The thrust chamber of the Space Shuttle Main Engine (SSME) consists of main 
combustion chamber with the injector plate assembly and the attached nozzle 
expansion section. Figure 1.1. presents the SSME combustion chamber and'the 
transonic flow nozzle. Figure 1.2 presents the SSME main injector assembly. 

The reactants: preburner combustion products, liquid oxygen and some liquid 

hydrogen, are injected to the thrust chamber through a main injector plate 
with 525 "concentric orifice" injectors (diffusive type burners) and 75 
"baffle elements" (premixed type burners). Propellants injected to the chamber 
are atomized, vaporized, mixed and combusted within the combustor volume. The 
reaction products are expanded through a subsonic/supersonic nozzle. The 
chemical kinetics of combustion reactions extends from the flame front to the 
expansion region. 

Proper design of a thrust chamber requires detailed knowledge of the flow 

pattern, temperature, concentration and pressure profiles as well as global 
parameters such as specific impulse, thrust, and wall heat flux. 

A mathematical model can provide such information in a cost effective manner 
and can aid in analyzing several designs and operational problems, for 
example: 

1. prevention of hot spots along the chamber wall and injector; 

2. prediction of pressure distribution and flow pattern (shock waves, 
boundary layer buildup and/or blowdown and/or creation of recirculation 
or stagnation zones); 

3. transient flow simulation (for start-up, shut-down and throttling 
situations); and 

4. influence of injection nonuniformities on the fluid flow, heat 
transfer and combustion within the chamber. 

Currently existing mathematical models for predicting the performance, propellant 
flow and combustion processes in rocket engines are identified in the "JANNAF 
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SSME POWERHEAD COMPONENT ARRANGEMENT 
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Rigorous Analytical Procedure" (Reference 1). In this procedure the 
performance of the thrust chamber is based on the application of the following 
two computer codes: 

TDK - "Two Dimensional Kinetics" (Reference 1) program for 

predicting the inviscid flow field using the method of 
characteristics with finite rate chemical kinetics, and 

BLIMP - Boundary Layer integral Matrix procedure (Reference 2) 

for predicting viscous effects in the boundary layer zone 
near the chamber walls. The program solves parabolic 
equations for momentum, and energy equations, with an 
algebraic turbulence model. It assumes chemical equilibrium 
of combustion reactions. 

The rigorous approach requires iterative calculations between these two programs 
in order to predict correct momentum and energy losses due to the interactions 
of inviscid flow and viscous boundary layer. 

Recent progress in Computational Fluid Dynamics (CFD) has lead to more advanced 
procedures and computer codes which are capable of simulating flows with above 
mentioned features. These codes solve Navier-Stokes equations by iterative solution 
algorithms. Since these equations represent fluid flow fields in general, solutions 
may be obtained for many conditions that exceed the existing TDK/BLIMP capability which 
considers only axisymmetric thrust chambers with cone or common bell shaped nozzles. 

Some future benefits from the advanced codes are: 

• Calculation of the subsonic, transonic and supersonic flow field 
with the same analytical concept for viscous and inviscid conditions. 

• Simulation of shock waves. 

t Flow separation at the wall, 
t Other than axisymmetric nozzles, 
t Dual throat thrust chambers. 

• Thrust chambers operating in wind tunnels or with ejector systems. 

• Flow over gaps or around protrusions in the nozzle wall. 

• Transient flow simulation for start-up, shut-down and throttling 
si tuations . 
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One such procedure has recently been applied to a three-dimensional, elliptic, 
two-phase flow problem to study the effects of rocket injector anomalies on 
wall heat transfer (Reference 3). Another example of new codes is CHAM's 
general purpose CFD code PHOENICS (Reference 4). PHOENICS has been applied 
to numerous diversified problems including several for NASA (References 5 to 8). 
Because of the flexibility of PHOENICS, and its current use at NASA MSFC for 
other projects, PHOENICS has been selected for the present study. 

1.2 Project Objectives 

The objective of this project was to demonstrate the suitability of PHOENICS 
(an existing general purpose multidimensional, elliptic, viscous flow and 
heat transfer computer code; Reference 9)., for predicting flow with combustion 
in SSME thrust chamber. 

Another important objective was to describe advantages and shortcomings of 
the proposed approach so as to assist NASA personnel in making comparisons 
between the CFD approach and the JANNAF procedure. Some of the pertinent 
questions are: 

1. How much time is necessary to prepare the input data to the code? 

2. What type of computer will execute the program efficiently? 

(Identify core requirement.) 

3. What is the computer program execution time? 

4. ' What is the accuracy of the analytical results? 

5. What is the status of the existing code and what problems pose 
significant difficulties? 

1.3 Outline of the Report 

The next section (Section 2) provides a brief .description .of.the .salj.ent . 
features of PHOENICS. Details of the specified SSME nozzle test case, selected 
computational grid distribution and boundary conditions are specified in 
Section 3. Implementation of this case in PHOENICS is described in detail in 
Section 4. 
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Computed results and their discussion are presented in Section 5. Specific 
answers to the above mentioned questions are provided in Section 6. 
Conclusions and Recommendations are provided in Section 7. All references 
are listed in Section 8. 
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Section 2 

DESCRIPTION OF THE PHOENICS CODE 


2.1 General Information on' PHOENICS 

PHOENICS is a general purpose CFD code capable of solving transport equations 
in arbitary geometrical configurations. The name PHOENICS is an acronym: it 

stands for Parabolic, Hyperbolic Or Elliptic Numerical -integration Code .Series. . 

The code consists of three modules: 

SATELLITE - for problem specifications; 

GROUND - for i ncorporation of new physical models; and 

EARTH - for solving the basic equations of fluid flow, heat and 

mass transfer. 

The SATELLITE program prepares data files, and acts in accordance with the 
information provided. Figure 2.1 presents the flow chart of the operation of 
the PHOENICS. It can be seen that there is only one-way interaction between 
SATELLITE and EARTH, and a continous interaction between EARTH and GROUND. For a 
more detailed description of EARTH, SATELLITE and GROUND the "PHOENICS Users 
Manual" (Reference 9) is recommended. 

Other salient features of the code are summarized below. 

• GEOMETRY 

The code provides great flexibility in geometry and grid selection including 

a) Body Fitted Coordinates (BFC) for representing complex domains with 
general nonorthogonal grid; and 

b) Porosity-Resistivity concept for representing blocked volumes, 
flow obstacles, perforated plates, etc. 

Both cartesian and cylindrical polar coordinates can be used. The available 
geometrical options offer versatility of Fin ite-Element Methods (FFM) and 
simplicity of Fdnite difference Methods (FDM). PHOENICS uses a Control 
Volume Method for discretization of the differential equations which, for 
regular geometries, is equivalent to FDM, and is briefly outlined in Section 2.3 
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It should be pointed out that the solution task in PHOENICS is similar to that 
in FDM as the coefficient matrices are positively defined and regularly 
diagonally dominant (in-FEM they are sparse and not always diagonally dominant). 
Both storage and computational effort are significantly greater in Finite 
Element Method (FEM). 

t PHYSICS 

I 

PHOENICS can solve up to 25 dependent variables. The specified dependent 
variables are: 

- velocity components for each of two phases; 

- enthalpies, temperatures and pressures for both phases; 

- two-turbulence quantities; 

- four concentrations; 

- the volume fractions of the interpenetrating media; 

- three equations for radiation fl ux components ; and 
few user specified variables. 

« NUMERICS 

Fully conservative and implicit formulations are used. As a result, 
there are no stability constraints, as commonly experienced in explicit 
and semi-explicit time marching methods. Equations are solved in 
"fully conservative" form ensuring both local (grid cell) and overall 
conservativion of all properties. This is difficult to accomplish in 
Finite Element Method or high order numerical methods. 

t PORTABILITY . 

The code is written in ANSI standard FORTRAN operationaT'on a large 
variety of computers; from mainframes as CRAY to mini computers as 
Perkin-Elmer , VAX, etc. Efficient storage management ensures that the 
computer space is only used for those variables which are solved for, 
or reserved by the user, at locations equal to the grid cell number. 

No change in COMMON or DIMENSION is required. The code dynamically 
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reserves the storage from one run to another. PHOENICS is user- 
oriented, well documented and verified for several problems. There 
are several "PHOENICS Demonstration Reports" (PDR's), from which 
potential users can benefit by adding more complex physics or geometry. 


2 . 2 Basic Equations 

PHOENICS solves the discretized versions of the basic differential equations which 
express the physical laws of "conservation" of mass, momentum and energy. The 
code has provisions of simulating both single and two-phase flows. In the two- 
phase flow mode, for a general conserved property, <j>, the transport equation is: 


aT 


+ div (rpV(£> - rf^grad 


<p ) = ifi + rS 


4 > 


( 2 . 1 ) 


where <j> stands for any of the dependent variables, r is the volume fraction of 
selected phases, p and T^ are density and exchange coefficient; V is the velocity 
vector, rfi is the interphase mass transfer rate (e.g. evaporation rate) and 
represents the source term. Equations for the various quantities differ primarly 
in the way in which the source terms and transport (exchange) coefficients are 
calculated. 


The above equation takes the special (simplified) form for the continuity equation 


(rp) + div (rpV) = r ft (2.2) 

The momentum equations are solved in each coordinate direction for each phase. 

For a single phase flow, in equation 2.1 and 2.2, r becomes unity and rfi becomes 
zero . 


A two equation turbulence model, is employed to calculate the turbulent 
exchange coefficients. 

Details of the k'vc model can be found in References 10 and 11. 
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Outline of Solution Procedure 


In Body Fitted Coordinates (BFC) the general finite difference equations are 
devised without any restrictions on the type of coordinate transformations, 
through the integration over a control volume. The resulting finite difference 
expressions are in conservative forms based on the staggered grid systems. 


The details of discretization in BFC are tedious and would require extensive 
vector algebra introduction and will not be derived here. However, most of the 
fundamental ideas stem from the procedure based on Cartesian coordinates in 
the original work by Patankar and Spalding (Reference 12, see also Reference 13). 

The resulting algebraic equations obtained from integration of differential 
equations over the control volume are of the form: 

a p4>p = d = E, W, N, S, H and L (2.3) 

where subscripts P, N, W, W etc. denote the grid point P and its neighbors N- 
north, E-east, etc. All coefficients are positive and the coefficient matrix 
is diagonally dominant i.e.: 

a p >_ £a d d = N, S, E, W, H and L (2.4) 

In the above formulae, a's represent the influences of convection and diffusion 
processes across the control volume faces. 

EARTH assembles all a's and efficiently solves the system of the algebraic 
equations (2.3) for all grid points and for all dependent variables. 

The solution procedure employed is based on the SIMPLE algorithm developed.) 
by Patankar and Spalding (Reference 12) in which the continuity equation is 
solved in terms of pressure corrections. The latest form of SIMPLE has been 
incorporated into PHOENICS. Another factor which influences the ability 
of PHOENICS to handle fine-grid problems, is the order of visitation which 
is selected when variables at cell-center are updated. 
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The order chosen involves what is called "repeated z-direction sweeps" through 
the integration domain (except for parabolic flows, for which only z-direction 
sweep is required). The whole set of cells is regarded as consisting of one-cell- 
thick "slabs", extending in the x and y directions, and piled one on top of the 
other in the z-direction. 


A single sweep therefore starts with attention being paid to the bottom (low-z) 
slab of cells. The finite-domain equations are solved for all the cells in 
this slab, the values of 4>'s obtained at the next higher slab being regarded 
as known. Attention then passes to the second slab, the ^-values there being 
adjusted by reference to those in the slabs both above and below. Then the 
next higher slab is attended to; and so on, until the adjustment sweep has 
been completed. For the solution of pressure corrections, two options viz: 
a) slabwise solution, and b) the whole field solution, are provided. 

Depending on the number of grid cells and problem nonlinearity up to a few 
hundred sweeps may be required before a converged solution is obtained. 
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Section 3 

DESCRIPTION OF THE TEST CASE AND NUMERICAL MODEL 

For demonstration of the computer code capability NASA MSFC specified the 
SSME geometry and operating conditions pertinent to the 100% power level. 

This section presents test case specifications, selected grid distribution, 
boundary conditions and assumptions of the mathematical model. 

This section presents a discussion of geometrical aspects of the calculation 
domain and assumptions of the mathematical model. 

3.1 Geometry of the SSME Thrust Chamber and Computational Grid 

Figure 3.1 presents details of the thrust chamber geometry with dimensions. 

The wall coordinates of the expansion part of the thrust chamber have been 
provided in Reference 2 and in data specified by NASA. 

A two-dimensional, axisymmetric , nonorthogonal body-fitted computational grid 
of NZ : NY = 41:20 control volumes has been selected for the computations. 
Figure 3.2 presents the computational grid. Figure 3.3 shows details of the 
grid distribution within the throat region. 

A nonuniform grid distribution has been used with finer grid spacings in the 
regions of steep gradients viz: in the throat region, and near the chamber 

wal 1 . 

The grid has been generated by employing the power law, which in the radial 
direction, is as follows: 

= 1 - (j|)“ j - 1. 2, . . . N - 1 (3.1) 

with y(N) = 1 and a = 1.5. 

Similar grid spacing variation has been obtained in the axial direction to 
concentrate the grid in the throat region. 
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Figure 3 . 1. Details of the SSME Thrust Chamber Geometry 




















GRID DISTRIBUTION NZ:NY (41:20) 



Figure 3.2 Grid Distribution of the SSME Thrust Chamber 


GRID DISTRIBUTION IN THE THROAT REGION 



Figure 3.3 Enlarged View of the Grid Distribution Within the Throat Region 
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An extra column of grid cells IZ - 41 has been placed downstream of the nozzle 
exit to accomodate fixed pressure exit boundary conditions. This practice will 
permit for more accurate calculation of the radial pressure gradients at the exit 
of the nozzle. It should be pointed out however that for better resolution near 
the exit a few extra slabs would be required. 


3.2 Physical Models and Assumptions 

The fluid flow, heat and mass transfer processes are described in terms of 
first principles for the entire computational domain. These are briefly 
discussed below. 


HYDRODYNAMICS 

The Navier-Stokes equation considered here in its full elliptic steady two- 
dimensional form is: 

V . p V = 0 ■ (3.2) 


V . (p V V 


?) = - vp 


(3.3) 


where c = -P ( 


9ui 

9xj 


9uj.\ , 2 
3xi' 3 


9ui 

3xj 


<5ij 


is the stress tensor. 


The momentum equations are transformed to general nonorthogonal form, similar 
to that of equation 2.3, and integrated within each control volume. 


The mixture density is calculated from the equation of state: 


p . M 


where 5 V 21 © 


.1 


j = H 2 , 0 2> H 2 0 


(3.4) 

(3.5) 


and where K is the gas constant, M is the mixture molecular weight, m. is the 

J 

mass fraction and M. is the molecular weight of specie j. 

0 


The effective viscosity in the momentum equations y is calculated as a sum of 
laminar, y and turbulent, u., viscosity. 

y = V* + Vt 
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TURBULENCE 


The Cebeci - Smith turbulence model (Reference 37) used in the past thrust chamber 
calculations is a "two-zone" model which accounts for the pressure gradients 
in the viscous sublayer and employs a Prandtl model in the outer boundary 
layer zone. For the purpose of this calculation, however, a more general, 
high Reynolds number k've. turbulence model of Launder and Spalding (Reference 10) 
has been used. The assumptions of Cebeci - Smith models for the viscous 
sublayer have been used to derive new wall functions used in model, which 
account for the pressure gradients. Details of the derivation and analysis 
are provided in Appendix B. 


The turbulent viscosity is calculated as: 
k 2 

U t ~ C d p £ (3./ 

where C^ = 0.09 is an empirical constant. 

HEAT TRANSFER 

For heat transfer calculations, the total energy (fr) is used as the dependent 
variable in energy conservation equation. 


For a premixed fuel and oxidizer composition, total energy is defined as the 
sum of enthalpy and kinetic energy: 




(3.8) 


where = nij /C pj dT + h f298j 
298 


(3.9) 


In the above relations summation is taken over all species, w is the local 
velocity, C^ ^ is the specific heat and is the enthalpy of formation. 


The specific heat C^ j is calculated as a function of temperature as follows: 

£ 

= a + a 2 T + a 3 T 2 + a 4 T 3 + a 5 T 4 (3.10) 

where a^ . . . a^ are constant coefficients given for each mixture component for 
two temperature ranges 300 - 1000K and 1000 - 5000K. 
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For an adiabatic flow, ft is constant throughout the calculation domain and 
therefore the energy equation does not have to be solved for. Temperatures 
can be calculated from algebraic relations (Equations 3.8 and 3.9). For a 
non-adiabatic flow, the wall heat transfer rate is determined from the Chilton- 
Col burn form of the Reynolds analogy: 

(3.11) 


q" = St 
^w 


P w 


- *w> 


where fr 

* p 

w 


P 

St 


w 


is the total enthalpy at the grid node in question; 

is the enthalpy corresponding to the prescribed wall temperature; 

is the mixture density; 

-2/3 

= Pr is the Stanton number related to the friction 

coefficient C^; 

T w 

= — ? ; and (3.12) 

ptr 

is the wall shear stress calculated from the wall functions (see 
Reference 10; see also Appendix B). 


The energy equation is solved in the following form: 

div (pVfr - Tgrad fr) = (3.13) 

where V = u eff /Pr is the transport coefficients and Sjv is the source term which, 
in general, should include the radiation energy transfer term and dissipation 
function <p. 

<f> - - V . (c ■ V) (3.14) 

which represents the rate of viscous forces energy transfer to thermal energy 
i.e. mechanical to interal energy transfer (for details, see Bird, Steward and 
Lightfoot, "Transport Phenomena", page 313). 

In this study, however, both radiation and viscous dissipation terms are 
neglected. 


Rigorous combustion calculations would require a full chemical kinetics model 
with transport equations for all species (H^ , H 2 , OH, Oj, 0 2 , H 2 0, H0 2 , . ..). 
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At high temperatures and pressures, the reaction rates are so fast in both directions 
that the chemical equilibrium is reached. In this case local composition can be 
determined from the prevailing pressure and temperature based on the Gibbs function 
minimization principle. 

In the present demonstration project, however, the simple approach of a one 
step reaction: 

H 2 + \ °2 * H 2° (3.15) 


has been applied with global reaction rate equal to slowest limiting chain 
propagation reaction viz: 

M + H + 0 2 ^ H 2 0 + M (Rl) 

This reaction, is of primary importance for fuel -rich flames in the main flame 
zone (Reference 14). Downstream of the flame front, within the throat and 
nozzle regions, recombination reactions become dominant. The reaction rate 
expression employed in present calculations is: 

R f = 2 . (H 2 ) (0 2 ) * 1.2 10 9 exp (-800/T) (3.16) 

3 

The terms in brackets are molar fractions with units in kg. moles/m giving 
rates in moles/sec. 


The fuel mass fraction m^ is calculated by solving the transport equation: 


div (Vm H - rgrad m^ ) = Rf 


(3.17) 


The calculations of the remaining stable species concentrations and m H g 
are obtained from stoichiometric relations. ^ ^ 


g 

= 8 . (m H - m„ ) ; and 


\-° ■ '"'h 2 ’ im h 2 


m H 2 0 1 ' m H2 " m 02 


(3.18) 

(3.19) 


D 

where m H2 is the m |H2 concentration at "fully burned stage" (in current case 

* 


when m QX =0.), defined as: 

B ° 1 0 
m H2 m H2 ' 8 m 02 


(3.20) 
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where superscript "0" indicates condition at the inlet (prior to combustion). 

Alternative, more adequate, combustion models for thrust chamber calculations 
are discussed in Section 5 as recommended for future study. 


3.3 Boundary Conditions 

For the two-dimensional computational domain, the boundary conditions should be 
specified at: inlet, exit, nozzle wall and symmetry plane, for all dependent 

variables, including: 

4> = w - axial velocity 
v - radial velocity 
k - kinetic energy of turbulence 
c - dissipation rate 
ft - total enthalpy 

m H2 " mass f ract i on 

In the present calculations the boundary conditions are specified as follows. 


3.3-1 Inlet 


At the inlet to the combustion chamber uniform profiles for all dependent 
variables have been assumed. These are specified as follows: 


total pressure 
radial velocity 
axial velocity 


P T = 2935.7 psi = 202.4 . 10 5 N/m 2 
v = 0 

w is calculated by the code based on the 
specified fixed inlet pressure 


enthalpy 

species composition 


ft = £ (C p T t + h°) 

ra, , m n , m u n are calculated from given 
u 2 m 2 u 

inlet mixture ratio f = 6.054851 and 

cal 


given inlet enthalpies h H = 1837.660 
- 3.84695 . 10 6 ^ 

h = - 2884. 385 3.77386 . 10 5 ^ 
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Details of the calculation for the inlet 
properties (m m m and T) 1[|let 
are provided in Appendix^A. 

k and e Initial values for k and e are estimated on the 

basis of inlet turbulence intensity of 10% 
and, length scale of 0.01 times the inlet 

chamber diameter. For the high-Reynolds - 

number flow in question, results are not 
expected to be sensitive to the inlet k and 
e values. 

3.3- 2 Symmetry Plane (Axis) 

At the symmetry plane the radial velocity v is set to zero while for all other 
dependent variables a zero normal gradient has been imposed. 

3.3- 3 Exit 

For hyperbolic or parabolic equations (solved by TDK and BLIMP code respectively) 
no exit boundary conditions are required. For fully elliptic calculations, 
however, a downstream boundary condition is necessary. A fixed exit pressure, 
uniform across the flow direction, has been imposed at the last slab (just 
downstream of exit plane). Although an approximate static pressure gradient 
for the SSME exit plane has been provided, it was decided not to use it. For 
the present project it would have increased the accuracy of the results in the 
exit plane domain. However, in a prediction mode for other thrust chambers 
this information is normally not available. Furthermore, the assumption of a 
constant pressure as a downstream boundary condition would yield results 
indicating, how far upstream this assumption would influence the calculated data. 
For the remaining dependent variables "locally parabolic" assumption at the exit 
plane has been, employed whj.ch. permits to calculate f-values . at the last .slab 
from upstream <j>-values. 
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Figure 3.4 Exit Boundary Condition Specification 


3.3-4 Solid Wall 

At the solid wall "no slip" boundary conditions are imposed. 

Special attention is required for calculating "near wall" shear stresses and 
heat fluxes where transport properties as well as dependent variables often 
vary steeply in the neighborhood of a wall. • For this purpose semi-analytical, 
solutions, called "Wall Functions" are used. 

Appendix B presents the basic principles of the wall functions for compressible 
flows with significant pressure gradients. 


3.4 Flow Field Initialization 

The starting conditions for steady state calculations are arbitrary. A "good" 
guess, however, can significantly improve the convergence rate. For the purpose 
of this study the initial flow fields have been generated from the isentropic 
flow relation for the nozzle (Reference 15). 


/A vZ _ , 2 , y-1 m 2, 

' 2 yfT M } 


M 


2 ) y-1 

Y+l 


(3.21) 
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For a given y = 1.3 and local area to throat area ratio one can calculate local 
Mach number M Q and then the local velocity, total temperature, total pressure 
and density. At the throat where A/A* = A / A t ^roat = 1 ^ ow is assumed t0 
be choked. A uniform total enthalpy equal to the inlet enthalpy and a uniform 


m 


H2 


mass fraction equal to the "fully burned" equilibrium value in the combustion 


chamber have been assumed throughout the domain. 


3.5 Under-Relaxation 

In the steady state calculations the sweep-to-sweep solution variation is 
controlled by the under-relaxation of dependent variables. The solution control 
is exercised by relaxing the pressure, velocity and density. Pressure and 
density are relaxed by the direct relation i.e. 


= 


n-1 


p* a + p 

p n = p"' 1 + p 1 a 


(l- % ) 


(3.22) 

(3.23) 


where: a - relaxation parameter; 

p* - currently calculated density; 

p n ~l ^ p n-l _ p rev j 0 us iteration values; 

p 1 - pressure correction. 


The axial and radial velocities are relaxed by adding an inertia term to the 

relevant finite difference transport equation: 

pVOL . * 
za.*.* + at. * 

♦p • W. 1 pVOL i - E, W, N, S, H and L, (3.24) 

1 A 4 * 


where <p stands fpr any jje pend enjt. vari_abl e; >t a. is ( the J.ink .cpefficient^ ^ 
(representing convectiv.e and diffusive. fluxes .through the. control.. volume faces); 
VOL is the volume of the grid cell; and Atp is the "false time" step used to 
control under-relaxation. Note that a smaller Atp value implies heavier 
under-relaxation. 
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In the present study the relaxation factors were modified during the calculations. 
For example the false time step Atp for velocity components were set to: 



• 01 • r Throat 
A * r Throat 
1000 r Throat 


for iteration 1 to 100; 
for iteration 100 - 300; 
for iteration 300 and up. 


Pressure relaxation has been gradually varied from a„ = 0.3 to a =0.6. 
- P P 
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Section 4 


INPUT DATA PREPARATION 

Input data preparation for the PHOENICS run requires user's input via the 
SATELLITE program and the GROUND subroutine. The input data are supplied in 
the form of FORTRAN statements. Both SATELLITE and GROUND are user oriented 
self-explana'tory subroutines with divisions into specif ic-purpose sections Land. . 
comment statements. In this chapter the description of the input data for the 
SSME NOZZLE test case is discussed. For general information about data 
preparation, the PHOENICS Instruction Manual (Reference 5) is recommended. 

Listings of the adapted SATELLITE and GROUND are provided in Appendix E. 

It is suggested that the reader who is not directly interested in the code 
execution procedure or in the studying of the coding implemented in SATELLITE 
and GROUND can skip the following two sections and continue onward from Section 5. 

4.1 Information Supplied Via Satellite 

The SATELLITE Program is divided into 43 Chapters of which Chapters 3 to 33 were 
used in the present demonstration case. Each chapter has the COMMENT heading 
and commented default variable specifications in <> brackets. A user may 
override the relevant variables. 


Grid and geometry are specified in Chapters 3 to 7 (see Figures 3.1 and 3.2). 

For the present two-dimensional case, in y-z (radial -axial coordinates) frame, 
XFRAC - nondimensional distance in ci rcumferential coordinate is set to 1 in 
Chapter 3. In Chapter 4 loop DO 400 is used to specify the nondimensional 
grid distribution in radial direction. NY=20 is a number of grids, POWER is 
the nonuniformity factor (P0WER=1 provides uniform grid) and THROAT = 0.13099 
is the throat region. The geometric shape of the nozzle is also specified in 
Chapter 5. In this case a "NOZZLE. DTA" file is read-in, containing (Z-Z^j iroa ^)/r^. 
and y wa11 /r T data. Current geometry data reflect the wall contour and is exactly 
the same as used in the BLIMP program.- NTAB = 393 is the total number of wall 
coordinates and NTABT = 31 is the throat coordinate location in the table. 

Array OZZLE (393, 2) contains these wall coordinates. In loops DO 500 and DO 501 
nondimensional ZFRAC-axial grid line coordinates are specified. Note that the 
zero location in the z-axis is at the inlet to the nozzle, not at the throat. 


4-1 



At the end of Chapter 5 an extra grid cell slab is added after the real nozzle 
exit and total number of grids is incremented by 1. 

In Chapter 6 the BFC = .TRUE, statement is set for nonorthogonal body fitted 
coordinate system. In DO 601 loop the arrays YN (NX, NZ, 1) and YS (NX, NZ, 1) 
are specified to north and south (wall and axis) y-coordinates of the nozzle. 

YN, used for nozzle wall coordinates, is interpolated from OZZLE array for 
each IZ, and YS is set to zero at each IZ since the south boundary is the axis of 
the nozzle. In double loop DO 605, DO 604 circumferential direction boundaries 
are specified as one hundred of a radian. In DO 605, DO '606 high ZH' and low ZL 
arrays are arranged such that ZL(IZ)=0 at the inlet and ZH(IZ)=ZFRAC (NZ)* 

THROAT is equal to the total nozzle length (including the extra grid cell). 

DEPENDENT VARIABLES - to be solved and stored are specified in Chapter 8. The 

indices of the SOL VAR and ST0VAR indicate selected variables i.e. PI - pressure, 

PP - pressure correction, VI - radial velocity resolute, W1 - axial velocity 

resolute, Cl - first concentration used for m u mass fraction, KE and EP for 

■ h 2 

k^e turbulence parameters, HI for enthalpy, H2 is used for temperature storage, 
U2, V2, W2 arrays will be used for storing the cartesian velocity components and 
storage 16 and 21 for nonorthogonal (contra-variant) velocity components and 
storage 23 will be used for saving continuity errors. Chapter 9 provides user 
specified titles to the variables. 

PROPERTIES - In Chapter 10, longest section of the satellite, physical properties 
of the test case are specified. SIGMA (24) and SIGMA (14) arrays contain the 
laminar and effective Prandtl numbers, EMULAM = 0.000102 is the laminar viscosity 
and IRH01 = -1 indicates that the density p will be specified in GROUND (also in 
Chapter 10). 

In the remaining part of Chapter 10 of SATELLITE, properties at the inlet to 
the nozzle are specified. The variable names are arbitrary and "local" in 
SATELLITE. FMIX = 6.054851 is the mixture fraction, SM0(1), SM0(2) and SM0(3) 
are used for initiating mass fractions of H 2 O, H 2 and O 2 , respectively. ENTHMIX 
is set to inlet mixture enthalpy as described in Appendix A. 

DO 1777 loop has been initially used to estimate inlet concentrations at 
TEMP = 300°K. The pertaining calculation procedure is also described in 
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Appendix A. Array SC contains molar concentrations of H 2 0, and 0 2 at inlet 

for T = 300°K; SC [kmol/kg], FMUB = 2 . SC ( 2 ) specifies the mass fraction of 

the "unburned" hydrogen and the next three statements override SC to represent 

a relative mixture composition in the fully burned stage. Note SC(3) = 0, 

i .e. m~2 = 0 as f < f .. . 

(J stoic 


CALL TEMPER passes the execution to a subroutine which calculates temperature 
TEMP for particular enthalpy HHH and gas composition SC. The variable TTEMPR 
is the initial "guess" of TEMP with RGAS = R = 8305 representing the gas 
constant and CPDR = C p /R the specific heat. All data are entered in SI units. 
CALL ENTHAL statement is used to calculate the mixture enthalpy ENTH = fr/R 
and specific heat CPDR = C p /R for given temperature TEMP and composition SC of 
three components. Both ENTHAL and TEMPER listings follow the listing of 
SATELLITE. 

WTMOL is the mean molecular weight of the "burned" composition. H2SAT = ft = 
RGAS. ENTH is used to transfer the total enthalpy ft to the GROUND where it could 
be used for "adiabatic" test case. If HI is solved for the H2SAT is not used. 


CALL MSOLV statement is used twice in the remaining part of this chapter. At 
first to calculate the inlet isentopic condition with given GA = y = 1.3 and 
AAT inlet to throat area ratio and secondly to calculate exit pressure. AM = 0 
and AM = 2 are the "guess" Mach numbers at the inlet and throat cross-sections 
which will be iteratively calculated in the MSOLVE subroutine. 

INITIAL FLOW FIELDS are specified in Chapter 13 via variable FI IN IT . Note 
that only the HI enthalpy is initialized as constant equal to H1SAT = H2SAT.TT0T 
(set just after the statement 6744). Remaining flow field variables PI, W1 , 

Cl and H2 (for temperature) are specified in a FLDDAT subroutine called at the 
end of SATELLITE, the listing of which is included in Appendix E •. F I IN IT = 
10101 indicates that particular variable will be specified in FLDDAT. Note 
that VI is not specified since the default value (l.E-10) is used. 

BOUNDARY CONDITIONS are set up in Chapters 14 and 15 by specifying the "REGIONS" 
at the inlet (Region 1) exit (Region 2) and nozzle wall (Region 3). No region 
specification is required for the axis because the symmetry "plane is a default 
region specification. 
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CALL PLACE (IREGN, TYPE, IXF, IXL , IYF, IYL , IZF, IZL) is used to specify the 
grid control cells in the region in question. 

CALL COVAL (IREGN, VARIABLE, COEFF, VALUE) is used to specify the boundary 
conditions for dependent variables. The boundary flux of the <J>-variable is 
calculated as: 

■» * CS m * C + ] (V 4 - *) (4.1) 

where <J> is the inner, adjacent grid node value to be calculated implicitly in 

PHOENICS, VALUE or v^ is the <j>-value at the boundary, C^ = COEFF is the 

diffusional flux coefficient of <b and S m is the mass source at the boundary, 

S m is calculated as: 
m 

S m = C m (V m - p) (4.2) 

m m m . 

where V m is either desired mass flow rate or external pressure depending on the 

type of the boundary condition, p is the pressure in the control .cell adjacent 

to the boundary (incorporated implicitly by PHOENICS) and C is the mass flux 

-10 m 

coefficient. C = 10 for fixed flow rate or a desired "resistance" related 
m 

value for the fixed pressure boundary condition. Detailed derivation and 
discussion of the boundary condition specification is provided in Appendix C. 
Note that ONLYMS = 0 indicates no diffusive link at the boundary. Note also 
that region 3 CALL COVAL statements are COMMENTed as the wall functions for 
all variables are calculated (in GROUND) by using a pressure gradient dependent 
formula . 

SOLUTION CONTROL AND RELAXATION parameters are specified in Chapters 26 to 34. 
For compressible flows, LOGIC (87) = .TRUE, is set in Chapter 26. FSWEEP = 401 
and LSWEEP = 450 indicate a restart run from results of 400 sweeps calculated 
prior to the current run. For the new run set FSWEEP = 1. 

In Chapter 29 the pressure relaxation factor FLXP and velocity under-relaxation 
parameter DTFALS are specified. This rather strong relaxation is reduced, in 
GROUND, after the first 50 sweeps. 

In Chapter 32 the logical arrays PRINT can be activated for whole field 
printout of desired variables. IZMON, IYMON are the grid coordinates of the 

4-4 



control volume values of which the (Wl, PI, T . ..) will be printed during the 
course of iterations (at each sweep). Variable NYPRIN = 2 will cause skipping 
of each second IY line from the output (NYPRIN = 1 can be used for full printout). 

RESTART run is specified in Chapter 42 where SAVEM = .TRUE', will activate result 
dump on the disk (file name TM1) which can be used for restart of graphical 
postprocessing. 

Finally, in the user independent section, calling subroutine FLODAT initial 
flow fields are prepared for all variables specified in Chapter 13. 

Additional subroutines attached to satellite are: 


SUBROUTINE TEMPER with input HSTAT - static enthalpy, TO - initial guess for 
temperature, SC (NSC) - species molar concentrations in [kmol/kg] and RGAS = 

8305 gas constant. The result of the calculations is T - temperature and CPDR = 
Cp/R. There is no user input required to this subroutine. 

TEMPER is called by SATELLITE and calls subroutine ENTHAL . 


SUBROUTINE ENTHAL - with input of: TEMP, and SC (NSC) returns mixture specific 

heat and enthalpy calculated as follows: 


^ = z 

R L \ 

+ + 

v 2 + 

v 3 + 

Z 5 T 

h 7 

, Z 2 T + 

Z 3 T2 + 

Z 4 T3 + 

Z 5 T 

RT Z 1 

+ T + 

3 + 

4 + 

5 


(4.3) 

(4.4) 


where h = 


lU 

f ,298 


/ T C dT 
298 p 


(4.5) 


Coefficients to Zg should be supplied via DATA statements in ENTHAL using 
the array ZA( 7, 2, NSC). Current data are valid for T = 300 - 5000°K for 
H 2 O, H^ and 0£ and are taken from JANNAF tables. 


SUBROUTINE MSOLV is called by SATELLITE and by FLDDAT to calculate isentropic 
flow properties. The input is GA = y = 1.3, KS = 1 for supersonic and KS = 0 
for subsonic flow, AAT - local to throat area ratio and AM - initial guess for 
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Mach number (e.g. Ma = 0 for subsonic and Ma = 2 for supersonic). The result 
of iterative calculations is the correct Mach number, and 


(4.6) 

(4.7) 

(4.8) 

(4.9) 

which are used in SATELLITE and FLDDAT to calculate local isentropic properties 

SUBROUTINE FLDDAT - is called from SATELLITE and shares with SATELLITE common 

block CMNBF1 (listing included after FLDDAT List). User dependent part starts 

in Chapter 2 where in two DO loops over IZ and MPH1 (variables) initial fields 

for PI, Wl, H2 = TEMP' and Cl - m u are set up. Note that within the loops CALL 

h 2 

MSOLV brings local cross-section isentropic parameters. Array PHI (IY, IZ) is 
used for initialization; 

4.2 User Supplied Coding Via GROUND 

The GROUND subroutine is accessed by EARTH at several stages of the solution. 

It has 16 chapters of which only the first 10 are used in the present test 
case. The COMMENT statements indicate when a particular chapter is being 
called by EARTH. User dependent section starts under the comment CXXX ... X 
USER SECTION 1 STARTS (see on page 1 of GROUND, Appendix E), where one must 
specify the array dimensions for CVAR, WAR, CM, VM and ZERO to be (NY, NX) 
i.e. for the current test case (20, 1). Any number of local variables and 
arrays, with proper dimensions, can be introduced in GROUND. It is advised 
to start FORTRAN names for arrays and real values with the letter G (GP - for 
pressure,, GW - for w-velocity, GTWALL .-.for wa.ll temperature, etc.) .This will . 
ensure no conflict between local variables and variables provided by GROUND 
station and EARTH COMMON blocks. 

Further downstream, a user can insert adequate DATA statement e.g. DATA RGAS/ 
8305.6/DATA ARFC , CRCF, ARCR CRCR/3. 7656E7 , 800, .../ where the second set of 


VRT = Ma ^ 

Y+l 

QRT = Ma l (y-) 2(y-1) 
TTT = 1 + ^ Ma 2 

T t 

PTP = (y^) 
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data defines Arrhenius reaction rates for forward (F) and reverse (R) reactions. 
Subsequently each chapter is briefly discussed. 

CHAPTER 0 - is called at the start of the run. Specific constants are set or 
calculated e.g. GPI = 4. arctg (1.) = tt . NXP1 = NX + 1, etc. Also file 23 
holding grid vertices (file name BFCXYZ) is opened. Finally subroutine WALDP 
is called to initialize some constants in it. 

CHAPTER 2 - in our test case will be called only once at the beginning of the 
run. Here the wall temperatures are calculated by calling SUBROUTINE TWALBC 
(discussed below) the input is NZ, GZNODE - nondimensional axial coordinates 
of the grid cell slabs, and ZTHRO - throat distance from the inlet. In DO 3120 
loop GZCELL (distance from the inlet) and GYWALL (local radius of the nozzle) 
are recovered from stored data of grid vertices on file BFCXYZ with Logical 
Unit 23. (These data were written by SATELLITE.) Finally wall inclination 
angles are calculated in DO 3140 loop. 

CHAPTER 3 - is used to gradually adjust the under-relaxation factors from sweep 
50 to sweep 150. This coding is test case depended and may not be optimal for 
the test case in question. It is assumed that more than 150 sweeps are necessary 
after which the relaxation is frozen. In the particular test case of the SSME 
nozzle approximately 400 sweeps are anticipated. Variables RLXP, RLXPZ, RLXPXY 
as well as Array DTFALS (NVAR) are EARTH variables available in GROUND. 

CHAPTER 5 - is used to specify source terms. As explained in Appendix C 
several physical processes can be simulated by providing appropriate "values 
and coefficients" of the linearized source term formula: 

SOURCE = [C m (V m - P) + C^] (V^ - <>) (4.10) 

In our test case C = V =0 and attention is focussed only 'on C.'and V, 

mm J <p cp 

(CVAR, WAR) coefficients with SOURCE = C^(V^ " 4>'). In the first part of 

Chapter 5, the reaction rate for m y is calculated as: 

h 2 

R = -A.e“ C/T . p 2 . m u m n (4.11) 

h 2 o 2 
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The linearization of above gives: 

R = R* + (m-,, - m *) = R* + R * " 0 


3m fu ' f u fu 


m,- * - m 


fu fu ,eq 


< m fu - m fu*> = 


•R* 


m fu " m fu ,eq 
Ae‘ 


(m fu,eq - "Vu 1 


“C/T p ni|_| nig 

m Ho* ' m H 0 „ ^ H 2,eq ^ 


(4.12) 


2 ,eq 


C. = CVAR 
<P 


V. = WAR 

<P 


where mu * - is the previous sweep fuel mass fraction; 
h 2 

m 0,, - oxygen mass fraction [m^ = f sto1c • (m fu * - !",„,«,)] ; 

A,C - Arrhenius rate constants set by the DATA statements as discussed 
above; 

p - mixture density (array RHO (IY,1)) 


Coding for the calculation of CVAR and WAR is executed in DO 508 loop. 


Note that CALL GET statements are used to "get" portions of the EARTH'S 
F.-array and "copy" it. to appropriate GROUND station array e.g.: 

CALL GET (C2 , CFU, NY, NX) 

accesses Cl = m^ u section of F-array and sets it into user's array CFU 
dimensioned (NY, NX) at the beginning of GROUND. 

After the CVAR and WAR are calculated they have to be placed into the 
appropriate position of the F-array. For this purpose the CALL ADD statement 
is used e.g. CALL ADD (Cl, 1, NX, 1, NY, VOLUME, CM, VM, CVAR, WAR, NY, NX) 
where Cl indicates the variable name, the following four integers define the 
grid cells range at the IZ-slab, VOLUME indicates that the source term is 
calculated per unit volume, CM = VM = 0 i.e. no mass addition is envisaged, 
CVAR and WAR are the calculated coefficient and value of the source and 
NY, NX are dimensions of the CM, VM, CVAR, WAR arrays. 
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In the remaining part of Chapter 5 wall functions are used to calculate the wall 
shear stress, wall heat flux and near wall k^e values. Between label 515 and 
517 local (IZ) near wall pressures, velocity and density are extracted from 
appropriate arrays GP, GRH. GPRL = SIGMA (24) is the laminar Prandtl number, 
set in SATELLITE. Four sections starting with labels 517, 520, 530 and 540 
call subroutine WALDP where wall functions are calculated. The input parameters 
are: IZED - current slab number, ISWP - .sweep number, NZ - total number of. slabs, 

1 - variable index 1-W, 2 - fr, 3 - k and 4 - e, GDYNY distance to the wall from 
the "near wall" grid node, EMULAM - laminar viscosity set in SATELLITE, GWFD1, 
GWFWl-w, two Prandtl numbers Pr^ , Pr^, GWFP1H - current slab and higher slab 
pressures, and finally, GDZ - AZ. The result calculated in WALDP is "VALUE" 
and "COEFFICIENT" set into CVAR (NY ,1 ) and WAR (NY ,1 ) and representing the 
near wall source term coefficients. Note that WAR (Wl) = 0 specifies a no 
slip condition at the wall and WAR (HI) = fr + Cp (T - T) is the wall enthalpy. 
The coefficients are wall shear stress and wall heat flux, respectively. Wall 
functions for k and e are used to calculated "fixed near wall" k and e. 

CHAPTER 7 - in the current test case calculations of dependent' variabl es 
including GATT, GFLXT - throat area and throat mass flow rate, GAEXT, GFLXE - 

exit area and exit mass flow rate are conducted for the last sweep only. 

In DO 772, DO 773 , and DO 774 loops axial components of the wall pressure 
forces are computed. Next thrust and specific impulse are calculated and 
printed. In the last part of Chapter 7 significant selected variables are 
printed. 

CHAPTER 10 - is used to calculate density and temperature. In DO 1010 loop 

species molar concentrations are determined from the known m u - mass fraction, 

h 2 2 2 

and molecular weight WTMOL. GEKIN is the kinetic energy equal to 0.5 (w + v ) 
and GHSTAT is the static enthalpy. The CALL TEMPER statement provides a new 
local temperature related to the* current* concentration' an’d static enthalpy. 

In the last three statements before 1010 label the density (from the equation 
of state), the temperature and specific heat are calculated. The listing and 
discussion of TEMPER subroutine is provided with the SATELLITE description. 
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CALL SET is used to set local arrays GRM and GTEMP into the F-array storage. 

In the second section of Chapter 10 selected axial and near wall pressures 
and temperatures are saved for later printout purposes. In the last sweep 
within the DO 1050 loop the Mach number is calculated for the entire flow field 
and set into the F-array via CALL SET (C2, . ..). Remember that the C2 - "second 
concentration" memory is being used to save and then print the complete Mach 
number field. 

SUBROUTINE TWALBC - is prepared to specify the temperature along the nozzle 
wall. User places ZTWQ and TQW arrays via the DATA statements. Within this 
subroutine a linear interpolation is employed to calculate the temperatures at 
the grid node locations from temperature data specified at locations along the 
wall. The User is advised to verify the resultant T^l array which is being 
printed from the DO 5555 loop at the end of this subroutine. 

SUBROUTINE WALDP - calculates VALUE and COEFFICIENT for the near wall source 
terms Wl, fr, k and e. The first section, up to statement no. 100 is accessed 
only once for initialization of prespecified constants. Between label 100 and 
120 local velocity.- WP, density RH0P and pressure gradient DPDZP are calculated. 
PPC is a constant in the p + ^ u + formula. 

RE is the Reynolds number calculated based on the control volume characteristic 
distance (cell size) near wall velocity and laminar viscosity. 

Further coding is divided into three parts: 

up to label 212 wall functions for w and fr are implemented; 
between label 300 and 400 wall functions for k - turbulent 
kinetic energy is coded; and 

after label 400 - the c wall functions are specified. 

The programmed analytical formulations for the wall functions are provided in 
Appendix B. 
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Section 5 


PRESENTATION AND DISCUSSION OF COMPUTATIONAL RESULTS 

The main objective of the present study is to demonstrate the capability of 
simulating reactive, compressible elliptic flows within the combustors and 
nozzles of a thrust chamber. This section presents the computed Navier-Stokes 
solution of the specified test case, already described in Section 3. 


5.1 Presentation of Results 

Figure 5.1 presents velocity vectors within the combustor and nozzle. Details 

of the velocity vectors within the throat region (IZ = 4 to IZ = 14) are shown 

in Figure 5.2 and velocity profiles at selected axial cross-sections are shown 

in Figure 5.3. Figure 5.4 shows Mach number contours with uniform contour levels 

distributed between Ma^ = 0.2 to Ma 2 Q = 4.0. The location of Ma = 1 contour 

is presented in Figure 5.5. Figure 5.6 presents contour lines of the absolute 

2 

value of static pressure p [N/m ]. The scale of pressure contour levels is 

4 7 2 

non-uniform and the range varies from p^ = 10 to p^g = 2.10 N/m . In Figure 
5.7 contours of absolute temperature T [°K] are shown. The range of temperature 
contour levels is uniform and varies from T = 1500 °K (Contour 1) to T = 4000°K 
(Contour 20). 

Figure 5.8 provides axial variations .of absolute temperature at three radial 
locations: 

T^ - wall temperature (input data); 

T^y - last grid cell near wall temperature (y/y w = 0.994); and 
T^ - temperature along the axis of the nozzle. 

Finally Figure 5.9 presents similar variations of static pressure in absolute 
values along the solid waif and along 'the axis'. ‘ *' *• 

Detailed numerical information of the flow field is provided in Tables 1 to 8. 

In Table 1 an output summary in SI units is presented where the headings have 
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Uocior scales 15625.0 m/3 

> 

Figure 5.1 Calculated Velocity Vectors Within the Combustor and Nozzle 


UELOCITV SECTORS IN THROAT REGION 
IZ - 4 TO 14,. 



Sector s cales 3125.0 m /3 

Figure 5.2 Details of Calculated Velocity Vectors Within the Throat Region 
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Figure 5.5 Location of Ma = 1 
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Rjinge of values is 6406.92 to 0 . 202I64E+08 [M/m 2 J 
Number of contour s>20 
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•PRESSURE CONTOURS P(N/M2) 



Figure 5.6 Calculated Pressure Contours 


TEMPERATURE CONTOURS TOO 
UNIFORM CONTOUR LEUEL DISTRIBUTION BETWEEN 

l-T-1500 AND 20-T-4000 
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Figure 5.8 Axial Variations of Absolute Temperature at the Wall, at the 
"Near Wall" Grid Node Along the Axis of the Nozzle 
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Table 3. Predicted Near-Wall and Axis Pressure Distributions (in Logarithmic 
Scale) 


IZ 

PW(PSIA) 

PA ( PS I A ) 

LOG ( PW ( PS I A ) ) 

LOG ( PA < PSIA ) ) 

1 

2 . 932 1 E+03 

2 . 93 1 0E+03 

3 . 4 67 2 E+00 

3. 4670E+00 


2 . 83 1 8 E+03 

2. 7628E+03 

■ 3 . 452 1 E+00 

3. 4413E+00 

3 

2 . 8 1 9 1 E+03 

2 . 6 9 92 E+03 

3 . 450 1 E+00 

3. 4312E+00 

4 

2 . 6786E+03 

2 .5711 E+03 

3 . 42 79 E+00 

3. 4 101 E+00 

5 

2. 3727E+03 

2 . 3580E+03 

3 . 3752E+00 

3 . 3725E+00 

6 

1 . 7 303 E+03 

2 . 0957E+03 

3 . 238 1 E+00 

3. 3213E+00 

7 

1 . 2807E+03 

1 . 8616E+03 

3. 1074 E+00 

3. 2 69 9 E+00 

8 

9. 590 IE +02 

1 . 6826E+03 

2. 9818E+00 

3 . 2260E+00 

9 

6. 1 4 69 E +02 

1 . 502 1 E+03 

2. 7887E+00 

3. 1 767E+00 

10 

3 . 2902 E +02 

1 . 2957E+03 

2. 5172E+00 

3. 1 125E+00 

1. 1 

2 . 3 7 5 2 E+0 2 

1 . 07 36 E+03 

2 . 3757E+00 

3 . 0308E+00 

12 

1 . 8 3 55 E +02 

8. 5263E+02 

2 . 263 3 E+00 

2. 9308E+00 

1 3 

1 . 4419E+02 

6 . 5 1 63E+02 

2. 1589E+00 

2 . 8 1 40E+00 

14 

1 . 1 460E+02 

4. 8330E+02 

2 . 0592 E+00 

2 . 6342E+00 

15 

9 . 22 1 9E+0 1 

3. 5 101 E+0 2 ' 

1 . 9648E+00 

2 . 5453E+00 

16 

7.51 86E+0 1 

2. 51 16E+02 

1 . 876 1 E+00 

2. 4000E+00 

17 

6 . 2327E+0 1 

1 . 7775E+02 

1 . 7947E+00 

2. 24 98 E+00 

18 

5. 1984E+01 

1 . 2508E+02 

1.71 59E+00 

2. 0972E+00 

19 

4 . 3898E+01 

8 . 7942E+0 1 

1 . 64 24 E+00 

1 . 9442E+00 

20 

3 . 75 1 8E+0 1 

6. 2038 E +01 

1 . 5742E+00 

1 . 792 7 E+00 

2 1 

8 . 2429E+0 1 

4. 4057E+01 

1.51 09 E+00 

1 . 6440E+00 

V 

2. 8322E+01 

3. 1 574E+0 1 

1 . 4521 E+00 

1 . 4 99 3 E+00 


2. 4999E+0 1 

2 . 2875E+0 1 

1 . 3979E+00 

1 . 3594E+00 

24 

2. 2273E+01 

1 . 6775E+0 1 

1 . 34 78 E+00 

1 . 2247E+00 

25 

2 . 00 1 3E+0 1 

1. 2461 E+0 1 

1 . 30 1 3E+00 

1 . 0956E+00 

26 

1.81 17E+01. 

9 . 3799E+00 

1 . 2581 E+00 

9. 7220E-01 

27 

1 . 6503E+0 1 

7. 1 573E+00 

1 . 2 1 76E+00 

8 . 5475E-0 1 

2 x 

1 . 5130E+01 

5 . 5 3 65 E+00 

1 . 1 7 98 E+00 

7. 4324E-01 

29 

1 . 3949E+0 1 
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9. 9981E-01 

7. 6209 E -01 

-8. 2532E-05 

41 

1 . 6494 E+00 

1 . 451 4 E+00 

2. 1 733E-01 

1.61 79E-0 1 
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Table 4. Predicted Wall Heat Transfer 
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Table 5. Computer Printout for Cartesian Component of Axial Velocity 
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Computer Printout for Cartesian Component of Radial Velocity 
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Table 7. Computer Printout for Absolute Temperature 
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Table 9. Computer Printout for Mach Number 
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the following meaning: 

ZG - dimensional axial distance (forward cell face) [m] 

ZND - nondimensional axial distance (Z - Zy)/ry[-] 

YN - local non-dimensional radius of the nozzle r/r T [-] 

ANG - wall inclination angle [rad] 

2 

PW - pressure near solid wall [N/m ] 

PAX - pressure near the axis [N/m 2 ] 

TW - wall temperature [°K] 

TNY - near wall temperature [°K] 

TAX - temperature at the axis [°K] 

Table 2 provides similar information in British Units p [psia], T [R°], where the 
heading (ZG:TH) = Z/ry represents nondimensional axial distance and YN - is the 
nozzle radius in meters. 

The wall and axis pressure distributions in the logarithmic scale (log P) are 
provided in Table 3. 

A summary of the wall heat transfer results is printed in Table 4 where 
2 

AREA (m ) is the wall area in contact with the control volume, HTCOEF is 

2 

the heat coefficient (W/m K) , TW-T(NY) represents the temperature difference 

0 ? 
in K between the wall and the "near wall grid node", QFLX (W/m ) is the 

heat flux, QDOT (W) is the local cell heat flow rate and QSUM ( 1 - I Z ) is 

the wall heat flow rate integrated along the wall from the injector face up 

to the IZ-th axial plane. 

Tables 5 to 9 are reproductions of the computer printout for the flow field 
variables including: 

Table 5, WCRT - cartesian component of axial velocity [m/s]-, 

Table 6, VCRT - cartesian component of radial velocity [m/s]; 

Table 7, TEMP - absolute temperature [°K]; 

Table 8, PI - absolute static pressure [N/m 2 ]; and 
Table 9, MACH - Mach number [-]. 

Calculated global performance parameters for the thrust chamber are: 

throat mass flow rate 

riiy = 456.72 kg/s = 1007 lb/s; 

mass flow rate at the exit from the nozzle 

m = 456.50 kg/s; 
ex 3 


5-17 



cross-sectional mass imbalance between the throat and exit 
m g = (456.72 - 456. 5)/456 .72 = 0.048%; 

thrust calculated as: 

T = Z P NY • A wall . sin (a wall ) = 19.353 . 10 5 N = 435,000 IbF 
summation is taken over thrust chamber wall and injector face. 

Figure 5.10 illustrates axial components of the pressure forces 
acting on the inner walls of the thrust chamber and used for 
calculating the thrust. 



Figure 5.10 Pressure Forces Acting on Rocket Chamber and Nozzle Walls 

Force acting on injector face was calculated based on the first node 
pressures rather than on the supplied data of inlet pressure. With 
the inlet pressure supplied the thrust is T - 440,000. 


specific impulse: 


I _ 43 5 , 000. ^2 o sec 

SP 1007 
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5.2 Discussion of the Computational Results 

The results obtained reveal all expected features of the nozzle flows. The 
velocity vector field presented in Figure 5.1 and 5.2 indicate flow acceleration 
in the expansion part of the nozzle. The viscous effects are visible from 
velocity profiles presented in Figure 5.3. As expected the velocity profile 
in the throat region has a maximum near the wall. While in further downstream 
cross-sections the max.imum velocity location is. shifted toward the axis. The 
viscous effects in the flow are visible in the downstream planes reflected in . 
the variation of the velocity profile extending up to the nozzle axis. 

The Mach number contours are shown in Figures 5.4 and 5.5 The sonic line 
(M =1) emanates from the geometrical throat wall point and reaches the 
nozzle axis at a further downstream location in a parabolic profile fashion. 

The shape and location of M = 1 line is typical for -the transonic nozzle 
flow field. 

Inspection of Ma = const contours downstream of throat region indicate the 
existence of a shock generated at the position of the second wall curvature 
tangency point (location where the radius of curvature meets nozzle bell shape). 
In this region Mach-isol ines exhibit a change of curvature in the vicinity of the 
wall. The "barrel type" shock wave is well visible in the region of fine grid 
just downstream of throat and extends almost up to the exit. 

The shape of the barrel shock is clearly resolved in the upstream part of the 
nozzle. The downstream resolution of the' shock is affected by the grid 
coarseness and by approximations in the exit boundary conditions specifications. 

Examination of temperature and pressure contours reveals that the largest 
pressure gradients exist within the throat and steepest temperature gradients 
downstream of the' 1 throat near the wall curvature tangency point. Note that 
temperature contour levels are presented with uniform intervals. 

In all contour-plot figures, the geometry outline extends up to slab 41, which 
is located downstream of the nozzle exit. The dotted line indicates the location 
of the exit. There are, however, indications of the influence of the uniform 
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exit pressure assumption on the inner flow field profiles. Contours of Ma , 

T and p indicate some "perturbations" near the nozzle exit. Therefore the 
results in the near-exit region should be interpreted with caution. 

Possible improvements in the specification of exit boundary conditions are 
discussed in Section 6. 

Interesting observations can also be drawn from Figures 5.8 and 5.9 presenting 
axial distributions of temperatures and pressures. The axial temperature 
profile T axi$ is typical for inviscid nozzle flows. The .near wall temperature 
T ny (in Figure 5.8) illustrates a complex phenomenon of expansion portraying 
the kinetic heating and convective heat transfer between the combustion gases 
and cooled nozzle wall . 

Distribution of the pressure profiles P ax -j s and P wa ]] s ^ ows that just after the 
downstream "tangency point" pressure near the wall becomes larger than that at 
the axis. Close to the exit the difference between the wall and axis pressure 
becomes: 

P ,, - P . - 39800 - 6900 = 32900 N/m 2 s 4.77 psi 

wall axis r 

This confirms our earlier conclusion, that for accurate calculation of "near 
exit" flow details, a better exit boundary condition is required. 

Finally, the analysis of the global parameters, such as total mass flow rates, 
thrust and specific impulse is summarized below. 

Cross-sectional mass flow rates are perfectly preserved 
throughout the chamber. 

Calculated, temperatures within the combustion chamber and 
the nozzle are overpredicted due to the "one-step" reaction 
model assumption used in present calculations. 

Finer grid in the regions, where significant axial pressure 
gradient exists is necessary for more accurate results. 

Predicted key parameters such as thrust and specific impulse 
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require an accurate pressure distribution in the contraction as well 
as in the expansion parts of the nozzle*. 

Further studies of the "boundary layer and main stream" interaction 
procedure is required with attention on "near wall" shear stresses 
and heat flux calculations. It is envisioned that direct and accurate 
calculation of the boundary layer structure is well within model 
capabilities and can help in understanding several phenomena of 
shock-boundary layer interaction. 

Application of a chemical "equilibrium model" and the complex 
kinetic effects are suggested for the further studies. The 
incorporation of finite rate chemistry would be especially 
interesting to study near wall reaction processes and interaction 
between chemistry and heat transfer. 


* Results of an additional run (results of which are not presented here) with 
four extra grid slabs in the contraction part of the nozzle revealed that the 
calculated specific- impul se I increased from 432.0 sec for 41 x 20 grid to 
437.9 sec for 45 x 20 grid. p 
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Section 6 

CFD MODEL ASSESSMENT 

The salient features and general advantages of the proposed CFD approach and 
computer code ( PHOEN I CS ) have already been described in Section 2. This 
section discusses the assessment considerations listed in the project 

objectives (Section 1.2, page 1-5). Each question is answered in specific 
reference to PHOEN ICS ' s performance for the problem considered. 

6. 1 Time Requirement for Input Data Preparation 

Section 4 has provided the details of the input data preparations. The time 
-requirement for input data preparation can vary depending upon the complexity 
of the test case and on the experience of the engineer. For experienced 
engineer modification of the existing SATELLITE and SUBROUTINE GROUND for 
the SSME thrust chamber test case may take: 

only few minutes for grid or boundary condition change or for 
additional printout arrangement, to 

few days to include new combustion model, new wall functions or 
other model improvements. 

To facilitate new users, CHAM conducts three-day training courses with hand-on 
workshop sessions. These courses are held once in every three months at CHAM. 
Specific courses are also held at client's site on request. 

6.2 Computer Requirement 

PHOENICS is a portable code. At present the code is operational at various 
main frames e.g. CRAY-1, CYBER 205, CDC-7600, IBM, UNIVAC, and 32-bit mini 
computers such as: VAX, PRIME, Perkin Elmer and Apollo. 

The code uses ANSI Standard FORTRAN and is not optimized or vectorized for 
any particular computer. For complex problems i.e. those with large number of 
control cells (e.g. _> 20,000), out-of-core storage is usually necessary on 
mini computers-. Due to extra 1/0 operations of the out-of-core storage use, 
this usually slows down the execution (up to 50%). The present work was done 
on a Perkin-Elmer, with 820 control cells and full in-core storage. 
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6.3 Execution Time Requirements 


The current problems required 400 sweeps (iterations) for the convergence of 
the iterative solution procedure. The convergence was judged by checking the 
mass residuals as well as the calculated flow variables. The sum of absolute 
values of residuals at all cells was reduced four orders of magnitude in 
400 sweeps. 

The sweep to sweep variation of flow variables were within 0.1% of their final 
values. 

Execution time of the base case was 90 * 110 minutes of P-E, 3205. Generally 
main frame computers are about 20 times faster. 

% 

6.4 Accuracy of the Results 

Accuracy of the results can be best determined by direct comparisons with the 
measurements. In the absence of measured data, the following considerations 
may be useful . 

1. There are no mass, energy or momentum imbalances in the calculations . 

This is due to the conservative formulation of the numerical scheme. 

2 . Results can be somewhat sensitive to the grid resolution. Generally, 
grid-sensitivity studies should be performed for each new problem 
This has been recommended for SSME Chamber problem, too. 

3. Results may also be sensitive to the physical models employed. Each 

of the physical models viz: for turbulence, combustion, wall shear 

stresses, and heat transfer have been described in detail. This should 
facilitate a fair assessment. Alternative models can be employed and 
assessed. 

6.5 Status of the Code 

The PH0ENICS code is in use at more that 100 research organizations over the 
world. Due to the extensive usage, the code is fairly well debugged. However, 
modifications for further improvements in the solution scheme as well as in the 


6-2 



user orientation of the code are in continuous progress. These are consolidated 
and released to the users periodically (i.e. once in a year or two ) by 
providing the upgraded code. 

The code has provisions for body fitted coordinates, which are essential for 
accuracy of nozzle flow calculations. Non-orthogonal grids can be generated 
fairly easily with the build-in "transfinite mapping" routine or can be 
read-in from user supplied grid preparation preprocessor code. 

The code has provisions for inclusion of new physical models by the user. 

This makes the code suitable for both routine engineering analyses as well as 
for research purposes. 
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Section 7 

CONCLUSIONS AND RECOMMENDATIONS 

Reported calculations have demonstrated that the basic capability exists of 
simulating reactive multidimensional flow in rocket engines for subsonic, 
transonic and supersonic compressible conditions by solving the Navier-Stokes 
equations for the entire nozzle geometry. Results of a single test case, 
presented in this present report, indicate good accuracy of the calculated 
data. All physical aspects of transonic flow with heat transfer in the nozzle 
have been predicted. Several improvements should be made in the mathematical 
model to make it a valuable tool for optimizing design and performance 
analysis of rocket engines. 

The proposed modifications are especially oriented to improve mathematical models 
of physical processes. With suggested improvements and after verification through 
Comparison with test results, sensitivity studies can be conducted leading to 
optimum rocket motor projections. 

Some of the improvements are discussed below. 

7 . 1 Improvements in Combustion Modeling 

The currently employed one-step reaction in the present study has been used 
mainly for convenience. Mimimum coding was required to calculate stable 
species concentrations (0^, and H^O). It is recommended that one of the 
following two techniques should be used for further improvement: 

a) chemical equilibrium model; or 

b) chemical kinetics of chain reactions. 

The first approach is the simpler one and requires no additional differential 

transport equation to be solved for. Local composition (O 2 , H^, 0, H, OH, H 2 , '^O, 

HO 2 etc. . .) can be obtained by solving a system of algebraic equations. General 
subroutines for the Gauss function minimization are readily available (Reference 
13) and can be incorporated into the calculation scheme. For the H 2 - O 2 
reaction mechanism a simple and fast procedure can be used without employing 
the Gauss function minimization. The task can be reduced to solving two 
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nonlinear algebraic equations for m^ and mg. The remaining species composition 
can be obtained either from the stoichiometric relation or from equilibrium 
constants. Appendix D provides detailed derivation of the above equations and 
outlines the appropriate solution procedure. 

The chemical kinetics calculations are recommended as a second step of the combustion 
simulation. The reaction mechanism between and 0^ is well established and 
rate constants are known. CHAM has extensive experience in using chemical 
kinetics modeling for the rocket plumes (References 19 and 20) o Relevant 
subroutines are readily available and only require careful incorporation 
into the PHOENICS 1 GROUND station and initial tests. 

7 . 2 Improvements in Wall Boundary Layer Treatment 

The "wall functions", applied in the present method, are based on experimental 
data for parallel flows in simple shear layers. The wall functions which account 
for the axial pressure gradient have not, to our knowledge, been extensively 
tested and require more attention during recommended verification studies. An 
alternative approach, also recommended for future investigations, has been 
recently presented by Launder (Reference 16) in which wall functions are 
replaced by local near-wall solution of the bounary layer equation. A separate 
grid structure is used between the near wall node and the solid wall. The 
resultant solution (for the velocity and for temperature) is integrated 
providing the wall shear stresses ( or heat flux) for the "basic grid" (macroscopic) 
calculations. A locally parabolic solution layer (PSL) is used to calculate boundary 
conditions for "macroscopic" equations. The proposed approach would combine the 
efficiency of implicit iterative schemes with the accuracy of fine grid time 
marching methods for shock-boundary layer calculations (see NASA Ames work in 
Reference 35 and 36). Figure 7.1 (Reference 16) schematically presents all three 
near wall treatments. 

a) currently used in PHOENICS "wall function" approach; 

b) fine grid solution of fully elliptic equations; and 

c) Elliptic - Parabolic Solution Layer (PSL) approach. 
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(a) (b) 



Figure 7.1 Illustration of Wall Boundary Layer Treatments 


7.3 Incorporation of the Thrust Chamber Cool ing System into the Computational 
Domai n 

In the present computational studies the computational domain extended up to 
the inner wall of the thrust chamber. The thermal boundary condition was 
provided in the form of specified wall temperatures T^-j. For the design or 
optimization studies, or even for transient analysis a priori specification of 
T t, is uncertain and inaccurate. The coolant temperature at the inlet to 
the cooling jacket can be specified more accurately. The coolant flow and its 
thermal coupling with the thrust chamber gases should be embodied into a general 
model with the computational grid extended up to the cooling jacket. 

'In future studies the coolant fluid flow and heat transfer as well as heat 
conduction within the thrust chamber walls should be included into the model. 

The necessary modi ficat-ion. woul d require specification of - 

- solid wall geometry (wall thickness, and shape) and thermal properties 
(thermal conduction and heat capacity); 

- cooling system geometry (inlet locations, jacket cross-section area, 
flow direction, etc), and; 

- cool ant properties . 
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CHAM has the necessary experience in simulating the conjugate heat transfer 
problems (References 33 and 34). 

7.4 Simulation with Two-Phase Flow Approaches 

Detailed information on the flame structure and combustion efficiency in liquid 
propelled rocket engines requires a two-phase flow and combustion formulation. 

The two-phases being considered are: 

- Liquid phase of cryogenic Hr, or 0^ streams, and 

- Vapor (gaseous) mixture of unburned reactants and combustion gases. 

There are t wo basic approaches to model the two-phase flow with combustion: 

Eulerian approach in which the gas and solid phases are treated 
as interpenetrat ing continuum, described by the conservation 
equations in the Eulerian frame, and 

Eulerian-Lagrangean (E-L) approach in which the liquid phase is 
described in a lagrangean formulation (coordinate system moving 
with the droplet) and the gas phase in Eulerian coordinates. 

The first technique has been developed by Spalding (20), Williams (21) and 
Westbrook (22) who modeled the spray through conservation equations utilizing 
a statistical distribution function in Eulerian frame. Similar approach has 
been employed by Swithenbank (23) for modeling a 3-dimensional gas turbine 
combustor. In all these applications a velocity slip between the gas and drops 
or particles has been neglected. An Eulerian approach for two-phase flows with 
velocity slip has been developed by Spalding (24) and used by Mai in et al (25) 
in two-phase coal combustion and by authors (26) of this document for spray 
combustion in diesel engines. 

The second (E-L) approach, introduced by Crowe (References 27 and 28) has been 
employed for several two-phase coupled fluid flows and combustion calculations 
(References 29, 30, 31 and 32). It has also been used for predicting the injection 
of fluid flow with combustion of rocket engines using liquid propellants (Reference 3) 
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The PHOENICS computer code has a built in option for the first approach (fully Eulerian) 
which can be easily used for the SSME combustion calculations. Specifications 
of appropriate boundary conditions and the evaporation rate formula are required. 

CHAM's past experience in diesel engine simulation (References 8 and 26) can 
be utilized effectively. 

7.5 Improvements in the Grid Generation- Procedure 

It is well known that the geometrical details of the nozzle have significant 
influence on the flow pattern in rocket engines. A Body Fitted Coordinate (BFC) 
system must be used for nozzle numerical studies. It is also known, however, 
that the grid selection can influence both the numerical convergence rate as 
well as the accuracy of the predicted results. Selection of the grid distribution 
for the entire calculation domain, in which the Navier-Stokes equations are 
solved, requires special attention. 

The currently employed nonorthogonal grid is generated by a simple method called 
"shearing transformation" in which radial Tines (z = const) coincide with cartesian y = 
const lines while axial lines are calculated at each z = const slab from 
nondimensional radial distance. 

An alternative practice would be to use an orthogonal grid, generated within 
the combustion chamber and nozzle. For tin's, however, a preprocessor grid 
generation package is required which would automatically generate the grid 
distribution. 

It is recommended that for the future studies of the nozzle flows a Grid Generation 
Package (GGP) should be developed to: 

a) quickly generate a grid in any two-dimensional domain; 

b) interact in any Navier-Stokes equation solution computer code; 

c) .> provide -resul ts n'n both computer interactive-session as well 

as "batch-job" mode; 

d) generate any coordinate system including orthogonal , nonorthogonal 
with easy specification of grid attraction in desired regions 
(walls, throat, shock waves, flame fronts, etc.); and 

e) be coupled with the N-S equation solver code in the "grid adaptive" 
mode. 
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CHAM has the necessary experience and personnel to quickly assemble, and 
deliver such a GPP code. Several phases of the BFC preprocessi ng have been 
already developed and are available as the GGP segments. 

7.5 Recommendations for Further Studies and Code Use by NASA MSFC Personnel 

PHOENICS courses and workshops, which are regularly organized at CHAM NA, will 
provide a good introduction to successful code utilization. In addition, it 
is recommended that the listing of SATELLITE and GROUND should be studied by 
the user. Most of the simple changes can be done by modifying SATELLITE. 

Familiarity with the code is best gained by performing calculations with 
different flow conditions. The next step can be to change the grid distribution 
in axial and radial directions. The following steps could involve calculation 
of fluid properties such as molecular viscosity as a function of temperature 
(modifications in GROUND ch. 11), calculation of laminar and turbulent Prandtl 
numbers, etc. 

Another action could involve the introduction of chemical equilibrium. A 
separate subroutine could be attached to GROUND and called from Chapter 10 
before the density is calculated. Appendix D provides a complete set of 
algebraic equations which can be solved by a Newton Method. 

For more complex modifications, consultation with CHAM is recommended to 
minimize the effort and hence cost and time. 
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APPENDIX A 


CALCULATION OF THE INLET CONDITIONS FOR 
SSME NOZZLE TEST CASE 



APPENDIX A 

CALCULATION OF INLET COMPOSITION 


For this study it is assumed that the fuel and oxygen are injected to the SSME 
combustion chamber in the liquid form. A single phase, homogeneous mixing 
assumption requires recalculation of the data provided for the liquid composition 
at temperatures below the boiling point T g to the equivalent composition at 
selected temperatures T. n above Tg. For the present study T >n = 300°K has 
been selected (one of the reasons is available C p % T dependence valid for 300 - 
5000°K temperature range). 

In reality SSME preburner combustion products are reacted with additional oxygen 
in a second combustion process in the main chamber. The assumption is made that 
the resulting reaction products are the same if oxygen and hydrogen are reacted 
in an equilibrium fashion with appropriate enthalpy values. 


The input data are: 
o Enthalpies 

ft u = -1837.660 ^4- = " 3. 84695. 10 6 J/kg 
H ^ mo i 


ft 


2884.385 ^- = - 3 .77386 • 1 0 6 J/kg 
mol 3 


Composition at the injection plane 

o 

m_ 


f = “7r~ = 6.054851 


m 


m 


H 0 1 + f 


= 0.141746 


H, 


m° + m° + m„ n = 1 
H 2 0 2 H 2° 


m 


o 2 1 + f 


m H 2 0 ■ 0 


= 0.858254 


Mixture enthalpy 


ft 0 = m° ft u + m° ft n = - 5.7768 -10 5 J/kg 
m H 2 H 2 0 2 0 2 3 
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/ \>Q Q 

Assuming an "adiabatic" process i.e. h = and T = 300 K the mixture 

composition can be calculated by a simple iterative procedure. By converting 

part of m° and m° im to m u n at constant the temperature of the mixture is 
o 2 n H^U m 

rising. During tne iteration process the g concentration was incremented to 
the level at which the mixture temperature reached T = 300°K. Calculated, in 
this way the inlet composition is: 

m H 0 = 0.0434175 
1- m H ? 0 

% = ~T~r~T - * 135592 

m n = m u . f = .82099 
°2 H 2 
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WALL FUNCTIONS FOR FLOWS 
WITH SIGNIFICANT PRESSURE GRADIENTS 



Appendix B 

WALL FUNCTIONS FOR FLOWS WITH SIGNIFICANT PRESSURE GRADIENTS 


Calculation of the transonic flow in the nozzles requires careful calculation 
of the wall boundary layer and main stream interaction. In the JANNAF standard 
procedure a few iterative computations between TDK - inviscid hydrodynamics 
and BLIMP - boundary layer flow codes are necessary to establish the correct 
boundary layer displacement thickness while matching imposed conditions at the 
potential wall . 

The simultaneous solution of the Navier Stokes equations for the entire flow domain 
implicity accounts for the coupling between the v-iscous boundary layer and the 
core flow. An accurate resolution of the boundary layer, however, is as important 
as in the JANNAF procedure. 

In the BLIMP program several grid nodes can be used across the boundary layer 
for accurate calculation of the shear stresses and heat fluxes. In the direct 
solution of the N-S equations the grid nodes have to be redistributed in the 
entire cross-section of the nozzle. To ensure ecomony and good accuracy in 
the direct neighborhood of the solid boundary a "wall function" concept is used. 

The wall functions permit to calculate the shear stresses and heat flux at the 
"near wall" grid node and the solid wall. A semi-analytical solution of the 
simplified N-S equations are used to devise the wall functions as discussed 
below: 


Implications of the Van Driest Formula 


The most commonly known Prandtl mixing length formula expresses the shear 
stresses x iin terms of the mixing length as 


X _ j>2 I 9u | 9u 
p 1 3y 1 9y 


or equivalent turbulent viscosity as 



9u 

9y 
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in which the mixing length £is assumed in the form of a ramp function (Figure 

B-l). 



where 6 is the displacement thickness, y - distance normal to the wall and x 
= 0.4 is the Von Karman constant. 

There have been numerous attempts to extend the mixing length formula £ - xy-. 
to include the viscous sublayer by modifying l. Van Driest proposed the following 

modi fication 

i • xy (i - e‘ y/A > B ' 3 

where A is a damping, wall constant calculated as 



and where v is the kinematic laminar viscosity. 


U e /x , is "friction velocity" and 

x w/p 


/\ + = 26 i s an empirical constant which may very with pressure 

gradient, transpiration, etc. 
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If we define a characteristic bounday layer distance y as 


+ 

y = 


- —L 

V V 
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(Note the similarity between y and Re), then the Van Driest formula can be 
expressed as 


Z = XY (1 - e 


+ . - + X 

-y /A ) 
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Various research workers have proposed further modifications of Van Driest 1 s 




exponent argument (arg = y /A ) viz: 


Patankar 

- Spalding 

+ 1/2 
arg = -y x 

Launder 

- Spalding 

arg = -y + x + /A 



+ * 

Cebeci 

- Smith 

arg = -y /A 


where A 


* = A + {?- (e 11,8u+ - 1) + e 11 * 8 ^} " 1/2 
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/ n + V dp + U 

T I = T/ V p = — d? ■ u =5- 

PU T 
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The Cebeci - Smith model extends the Van Driest formula to flows with significant 
pressure gradients. 


Van Driest formula can also be used to calculate an effective viscosity 


P eff = v + £ 2 ||y| = v + x 2 y 2 U - e 


-y + /A + ) 2 ,3u, 

3y 
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or in nondimensional form 


v 


+ . 1 t «_ |M| . 1 + x!i! (1 . e -y + /A + ) 2 ,3u, 

V 1 3y 1 V u hy' 
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Figure B-2 presents the resulting distribution v + ^ y + (see Reference 10). 
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Figure B-2 


For y + ^.11.6 a simple relation can be established 

+ + + n 1 9 

v -* v = xy D-i £ 

which will be used as a basis of wall funcitons 
Wall Functions 

Finite difference schemes used in fluid flow calculations approximate the 
velocity gradients simply by a linear difference 


du „ u j+l ~ u ,j 
dy Ay 
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Near the solid boundary where large velocity gradients exist the shear stress 
(expressed in terms of the velocity gradient) 


-t 


w 



u 
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is obtained by solving simplified N-S equations viz: 



B- 15 


d / du_\ _ T w d£ 

dy ^ p eff dy' dy dx 


which after integration yields: 


du 




p eff dy T w + y dx 
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Introducing nondimensional quantities: 


V 




u = u/u. 


+ _ V d£ 

p 3d7 

PU T ' 

+ + veff 

y = v = 


= 1 + -£ 
V V 
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u = /t /p" 

T WV K 


we can write equation (B-15) as 


■f + + 

du _ 1 + p y 
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, T T 

dy v 

First Case P* = 0 . 

If axial pressure gradient is negligible P + = 0, then integration of equation 
(B-18) with equation (B-12) i.e. v + = xy + gives: 


/^ T dy + = fKrr. dy + .=. I £n (Ey + ). * . 8-19 

0 dy 0 X y + X 

where E is the integration constant E = 9 calculated for u + = y + with y + < 11.6 
from 

+ + X . 4* V 

u = y = - fcn (Ey ) 

A 
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The expression below is the well known log-law wall function formula 


u + ' = ^ £n (Ey + ) 
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Finally the shear stress is calculated from the relationship 


, 3u + U wall - u 

T w /P = v eff 9y = v v Ay 
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where v + is obtained from eq. (B-18) with P + = 0 


, + = it, s 
. + + 
du u 
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Second Case P > 0: 


For flows with significant pressure gradients P > 0 an analytical solution of 
eq. (B-18) viz: 


, + - + + 
du _ 1 + P y 

, + + 
dy v 
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has been obtained by Launder and Spalding (Reference 10) in the following form 


u + = \ (2{/l+ Py - 1} + ln{ 

A 


4Ey 


+ + 2 A + P'y 


2 + P y + 
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where E is a function of P and can be calculated from the above formula with 
u + = y + = 11.6 


The shear stress is calculated similarly as in the first case. 


Implementation of the Wall Functions 

■j- 4" 

The previous equation for u is nonlinearly dependent on y and p . Therefore 
we will use definitions of y + , p + and u + expressed in terms of one common 
variable s - called the shear rate coefficient 


s 


T w _ T w 1 


2 2 
pu p u 


2 

U T 


/ +\ 2 

(u ) 
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yu 1 _ Re 

+ + 
v u u 

+ \3 


Re /s 


^-£- z (JL ±- c • <“ + ) 3 

0 u 


n + + r, i +s3 Re r n / +\2 c Re 

P y = C(u ) — = C • Re (u ) = — — 

u 


Before inserting them to u notice that: 


(1 + /I + py ,) 2 = 2 + py + 2 J\ + p y 


Lets define, for convenience 


, a = /l~ + p + 


p’y + then u + is written as: 


u + (2a - 2 + An (4Ey + ) - An (1 + a) 2 ) = 


= | An (Ey + ) 


- £ (2 - An 2) + ~ (1 + a - In (1 + a)) 

X. A 


standard wall 
function formula 




additional terms due to 
pressure gradient effects. 
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APPENDIX C 

THE IMPLEMENTATION OF BOUNDARY CONDITIONS 
IN PHOENICS (EXTRACTED FROM PHOENICS USER CLUB 
(PUC) NEWSLETTER #3). 



THE PHOENICS IMPLEMENTATION OF BOUNDARY CONDITIONS 


1 ) introduction 

Many people using PHOENICS have trouble with a boundary 
condition or source description when setting up their SATELLITE. This 
section will be concerned with the explanation of the conceots behind 
the PHOENICS boundary conditions. How the correct COEFFICIENT 
and VALUE are determined will be shown: and examples of frequently 
encountered boundary conditions will be given. It will be assumed that 
users are familiar with the subroutines PLACE and COVAL to specify the 
locatfon and the appropriate coefficient and value for the 
boundary/source. 

2) The General Principle 

All boundary conditions are Inserted in PHOENICS as sources or 
sinks of one or more variables. These are represented by linear 
expressions of the form: 

s* ■ it* ♦ CC S(* ]]} * (V* - *F) ( 1 ) 

where the mass source s m is given by: 
s. . C. (V. - Pp) 

The symbol tl )1 refers to the use of the upwind practice, the subscript 
P refers to the centre of the finite-domain cell in question, and C$. 
V$ are referred to as the 'coefficient' and the 'value' for variable <t. 
Expression (1) can be used either 'per unit area or volume' or as the 
total for the cell . In question. Equation (1) is easy to understand, if 
we recall the general form of the PHOENICS finite-domain equations, 
eg: 

*P *P * 1 *1 -*1 ♦ (Su - S P ®p) (2) 

Equation (2) was derived by considering the contributions to the 
balances over the control volume P. of the neighbouring control cells. 



For any Internal link. In a given control volume P. say the link 
between point P and Its East neighbour, we write the term: 

A E (*E 0 P) (3) 

where the coefficient Ag represents effects of convection and diffusion 
added together. Including the links of P with all Its neighbours, and 
collecting terms, we arrive at the general equation (2). 

In general, for Internal points, the interactions between any two points 
are expressed by adding the term: 

*1 Ink (^n« Ighbour “ ^p) f 4 ) 

-"•r« 8 i ink ■ 0 - CC"] 3 (S) 

D stands tor diffusion and Urn]) for upwinded convection. 

PHOENICS treats boundaries In the same manner as any other points 
In the field Therefore, the effect of a boundary on an adjacent cell 
is expressed by adding to the finite-domain equation a term: 

a 1 ink (*b * *P) < 4 > 

where <j>b Is the boundary 'value' of <t>. Again: 

a 1 ink * 0 ♦ CC* J J 

Therefore equation (6) becomes: 

(0 ♦ CC-]]) (* b * *p) ( 7 ) 


r i 




which is identical to equation (1) if: 



Again Urn]) symbolises that, in accordance with the 'upwind' principle, 
only inflow is of influence. Therefore, for the general variable <t>. the 
‘coefficient* is simply the diffusive transport of o between the boundary 
itself and the centre of the adjacent cell. 

The convective transport at the boundary [lm]} requires special 
treatment, because PHOENICS does not store the boundary values 
themselves. Define: 


i - C. (P b - P P ) W 

which simply means that the flow rate is the pressure difference 

(between an externa) pressure and the pressure at P) divided by the 
flow resistance. Thus C m has the physical significance of the 
reciprocal of. a flow resistance between the Internal grid node P, 
where the pressure Is pp. and the external pressure p^. which we 

shall denote as V m (for value for mass) for the purpose of retaining 
the general concept of 'coefficient and value'. 

3) The Flnlte-Oomaln Equation for Boundary Cells 

The finite-domain equation for boundary cells is the same as 

equation (2). where one of the neighbours is missing (the one 
corresponding to the boundary cell face) and where the linearised 

source S,j - Sp qp contains the additional source/sink due to boundary 
conditions. Thus, the equation becomes: 

( A P ♦ S P ♦ c $)*p “ A i p p ♦ S u * % C’) 

• i c ap t boundary 

The pressure-correction equation becomes: 

(*P * c.) Pp * C A. p' . C m V B * E ' ( 10 ) 

•leapt boundary 

where E is the mass error. 


4) The Treatment of Various Boundary Conditions According to the 
Above Formulation 

( l) Special Cases 

( a) Fixed Source of Mass. Regardless of Pressure 
Set: 

C. - 1 .£ - 10, V. • S. . 1/C. * 

where S m Is the desired source of mass. This Is a direct 
consequence of equation (8) since. 


* See note in Section 7 

» - % m * C. (V. - p p) » 1.E-10 . (S,, . 1 . £ 10 - Pp) 
* S* - 1.E-10 Pp = $„ 


(b) Fixed External Pressure Regardless of Mass Flow: 


Set: 

C. - i.Eio, v* * p- 

where p* Is the desired value of pressure. This case will be 
discussed further in the next section. 

(c) Fixed Value of Any Dependent Variable. 

Set: 

Cq • 1.E10, V 0 • 

where Is the desired value of <t>. This is a direct consequence of 
equation (9) since: 

C A i <t> 1 ♦ Su ♦ 1 . E 1 0 • o* 

* 4p 7 Sp ♦ l.lio 

The 1.E10 terms are much larger than the other terms and therefore: 


<fcp 


1 .E 10 « »* 
1 . £ 1 0 
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( II) Treatment of Various Practical Boundary Conditions 
( a) Linear Source (For Example Walt Heat Transfer) 


Sh * ^waii = a (Tya11 ' 7p} = a/cp {cp 7 *aii " h p) 
where a = heat-transfer coefficient. Therefore: 
th 3 a/c P ; v h 3 c P T ua)] 


< b) Non-Linear Source (For Example Internal Resistance) 

♦ Distributed resistance: 

S w * K dr pw ? A. 5 z ( 1 e dp/dz * K dr pw 2 ) 

Th#n: C w * K dr pw A Sz , V w * 0 

t Baffle: 

S y » R D pw 2 A ( i* pp • K b pw 2 ) 

Than: C y ■ K b pw A, V w * 0 

# Pressure-drop per row of rods In cross flow 

PkP> coll * K 1^2 pw 2 ti , where N is the nu«ber of rows 
Sw, total » * 1/2 pw 2 A N 
Then : C y • K 1/2 pw A N , V y * 0 

5) The Inflow Boundary 

Fixed mass flow (m); applicable In both single and two phase 
flows) . 

Set: 

c. ■ i.E-io, v p • m 

Also set V$‘s for all p's. this is because a Mow rate of m carries Into 
the domain the boundary values of alt P's; for example It carries In 
w-momentum equal to mw| n . energy. mhj n . etc. 

In circumstances where convective transport across the boundary Is 
dominant. C$ would be set to zero. This Is usually the case. If- It Is 
desired to Include diffusive, as well as convective, transport. Cp would 
be set to the appropriate finite value r/s. where r Is the diffusive 
transport coefficient and 5 the distance between the boundary and the 
first Internal grid node. 

6) The Outflow Boundary 

Usually the specification of a fixed-pressure boundary condition Is 
adequate. This is achieved by setting: 

C m - "large", V B - required p 9Xt 

This boundary condition setting has the meaning that the outflow Is 
discharging in a reservoir of uniform pressure p 0 xt and ,eads t0 P P 5 
(at the last grid node before the outflow) such that the expression 

m • Z m (p p - P #x t) 

provides the correct flow-rate (satisfying overall continuity) Users 
familiar with recent finite-element techniques would recognize this 
practice as the 'penalty method’. 

Caution must be exercised in choosing the 'large* value for C m - 
particularly when p Qxt is not equal to zero; for it is quite possible. If 

C m is too large, for the value of < p - P Qx t> < or the r0dUired outflow 

(or inflow) rate to be so small compared with the value of P 0X t as to 
be of the order of. or less than, the precision of the computer 

storage. This can result in large errors in the evaluation of ( p 
Pext > and hence of the flow rate, which can adversely affect the 

solution. 
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In general.' for the case of a fixed external pressure. C m would be 
prescribed lo give the required dependence of the flow rate on 
pressure difference. Thus, if C m were set to 1.0. then a pressure 
difference <p Q *t - p) of 1 Nm~ 2 would cause an inflow rate of 1.0 
kgs -1 or kgs~ x m~ 2 depending on the setting of type) . A good rule is 
to set C m to be of the order of: 

1000 * (expected mass-flow rate / p Qxt ) . * 


To specify a known resistance or pressure drop at outlet set: 


C. ' 


Ap 


V. 


P«xt 


The rule is slightly different for two phase flows. For an outflow (le 
Cm. I (v m.l ~ P J < 0) the outflow is multiplied by Rj and pj to 

become: 


Pi ^ 1 , i (V>, j • p)j 

where pj and Bj are the local values of density and volume traction of 
the phase. This ensures that, if the values of C m \ and V m j are the 
same tor both phases, the outflows of each phase from a cell are In 
proportion to the mass concentrations of the phases in the cell, as is 
physically reasonable. 

7) Note 

When using the subroutine COVAL several FORTRAN variables have 
been provided to simplify the input. These include: 

F IX FLU (- 1.E-10) 

F1XVAL (- 1.E10) 

OHLYHS (-0.0) 

If FIXFLU or 1 . E- 10 is used in the call to COVAL the multlpltcaton 
of the value by 1.E10 is done by COVAL and therefore should not be 
done by the user. 


* This formula is very approximate and the coefficient found may 
need further adjustment - by up to *,2 orders of magnitude to give the 
correct flow behaviour. 
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APPENDIX D 


THE EQUILIBRIUM CHEMISTRY MODEL FOR H £ - 0 2 SYSTEM • 

(Based on N.C. Markatos, D.B. Spalding and D.G. Tatchell 
"Combustion of Hydrogen Injected into a Supersonic Airstream" 

NASA CR - 2802, 1977). 



Appendix D 

THE EQUILIBRIUM CHEMISTRY FOR H 2 - C> 2 SYSTEM 

The equilibrium chemistry model used to predict the properties of hydrogen- 
oxygen flame can be expressed by the four following reactions. 


H + H •* 1 H 
k 2 

0 + 0 t 0 o 


D-l 

D-2 


H + OH i ° H 2 0 


0 + H j OH 


0-3 

0-4 


There are six species involved in these reactions and only two atomic elements 
viz H and 0. Mas balance relations provide two basic equations in the’ form of 
total atomic mass fractions for 0 and H: 


p _ m , A w o w 0 

F ~ m O m O W m H 0 + W — m nw 

U 2 u W H„0 ll 2 OH 0H 


w. 


Hr 


w 


H 


f *. m u + m H + rj — m „ n + w 

H 2 H W U 2 0 H 2° W OH 0H 


m. 


D-5 


D-6 


where F is the total mass fraction of oxygen in any form and f is the total 
mass fraction of hydrogen in any form. In H 2 - 0 2 flames the following 
relation holds: 

f + F = 1 D-7 

Therefore if F is known, f can be determined. 
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From thermodynamic considerations the 
Kp, for the reaction: 

K 

P 

aA + bB ^ cC is defined by, 
c 


K = 
P 


X C 


c-a-b 


a b 
X A X B 


equilibrium constant 


( 0 - 8 ) 


where x stands for concentration and the pressure p is in 
atmospheres. For each of the four reactions ( D-l ) to (D-4 ) in 
the present model c-a-b = -1. Expressing the concentrations 
in terms of mass fractions by noting that, 

Wi 

m i W x i ( D-9 ) 


where W is the molecular weight of the mixture we get: 


W 


K 


= K nW 


pP’ 


C 


Wo 
A B 




m 


b 

B 


(D-10) 


Thus the equilibrium equations for the reactions (D-l) to 
(D-4 ) can be written: 


I 





m 


0 


f 



"V 

m H m OH 


( D-ll ) 


( D-l 2 ) 


( D - 1 3 ) 


t 



m 


OH 


tT'H 


(D-14 ) 


The condition of equilibrium is expressed using four 
equilibrium constants for the four chemical reactions. If 
thermodynamic equilibrium prevails ( the K's take values 
which depend upon temperature alone. In the above set of 
equations there are seven unknowns (m^ ( , m^, m QH 1 
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and F; f being given) and seven equations (0-5, D-6, D-7 

0-11, D-12, D-13, D- 1 4 ) The problem is therefore soluble. 

The remaining discussion defines the solution procedure. 
Derivation of Solution Procedure 


The present procedure is based on the reduction of the 
number of variables under consideration. Two equations 
are derived as follows: 

Equations (0-11 to D-14) ,are solved to give the relationships: 






m. 






» 






r 



(D-15) 


(D-16 ) 


( D - 1 7 ) 


(D-18) 


It is convenient to define the following parameters which depend upon 
temperature alone and therefore, can be tabulated right at the begining. 



(0-19) 


C 


K 3 K 4 


» / f 
K 1 K 2 



I 



(D-21-) 


B 


( D- 20 ) 


D 


(D-22 ) 



Using equations (D-15 to D-1'8 ) to eliminate ni Q , m R Q> nr^ and 

2 

rriQ R from equations (D-5 ) and 0-6 ) gives: 


F = m + /m 

u 2 U 2 


<s + 1 c % + if B 


(0-23) 


f = m 


1 » 


H 2 (1 + 9 + ^ 


iA * jV C 


(D-24 ) 


where the definitions given by equations (0-13 to D— 22 ) have been used. 


These equations have the form of a quadratic equation: 
-2 

au + bu + c = 0 

which yield the solution: 


u 



4ac 


(D-25) 


Note that a and b are always positive and c is always 
negative. Since u > 0, the physically meaningful root is 
the one with the positive sign. This particular form of 
quadratic expression is chosen since it does not require 
subtraction and gives greater precision. Using equation 
(D-25 ) to express the solution of equations (0-23 ) and (0-24 ) 
gives : 


m 


0„ = 


B + 1 c n ii 2 + if c ^ 


" + f c m n 2 * if ° 


2 

+ F 


Equations ( D - 26 )and( D-27)contain the unknowns and m u 

°2 H 2 


which can be easily determined, now< 


(0-26) 


v> 


r i 


r 

r i ~ 

2 

5 * 17 D 

+ J 

1 

A + T7 D /m 0 2 


2 J 



L 2 J 


o. 


(0-27) 
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APPENDIX E 


LISTINGS OF SATELLITE AND GROUND 
FOR SSME NOZZLE TEST CASE 

Listings of SATELLITE includes: 


a) 

Main SATELLITE 

pages 

El 

to 

E7 

b) 

Subroutine TEMPER 

page 

E8 



c) 

Subroutine ENTHAL 

page 

E8 



d) 

Subroutine MSOLV 

page 

E9 



e) 

Subroutine FLDDAT 

pages 

E10 

to 

Ell 

f) 

Common Block 

page 

Ell 



Listing of GROUND includes: 





a) 

Subroutine GROUND 

pages 

E12 

to 

E20 

b) 

Subroutine WALDP 

pages 

E21 

to 

E22 



i":- — - GROUP 3 . X -DIRECT I ON 8 
C N X < 1 > , XU LAST-;' J. , 0> , XFRAC ( 1 -30 ) 

XFRAC < .1 ) - 1 . 

C GROUP 4, Y-DIRECTION : 

C MV< 1. > , YVI.AST < 1 .. 0> . YFRAC ( 1 -30 > , R I NNER . SMALFA 

NY- 20 
PO'-IRR- 1 . 5 
T HE 0 AT—. 1 3 0 8 8 
YF RAO. ( NY ) ~ t . 

DO 400 I '/= 1 , NY— .1 

400 YFRAC < NY- I Y )«1 . -( FLOAT ( 7 V ) / FLOAT < NY ) ) ** POWER 
WR I TE ( 6 , 284 7 > ( YFRAC < I Y ) , I Y« 1 , NY ) 

2847 FORMAT ( " YFRAC” / ( 1 P5E 1 1. . 3 ) ) 

C GROUP 5. 7. -DIRECT I ON s 

C ' NZ < 1 > , Z WLAST-C 1 . 0> , ZFRAC ( 1 -30 > 

C— M7.MAX - 47 

NZ-40 
MZT-6 
N FAB— 3'? 3 
MTABT--3 1 

C- -READ U'! NOZZLE GF.OI-IETRY DATA FROM FILE NOZZLE. UTA - 

OPEN ( 25 , F I L E- "NOZZLE . OTA , ACCESS* " SEQUENT I AL , STATUS- " OL D " 
i. , FORM-- -• FORMA ('TED " , ERR-99? ) 

READ ( 23 v 4 ) U l l LE 

ZFRAC ( NZ ) =0Z ZLE ( NT A 8 , 1 ) -07. 7.LE (1.1) 

Z F R A C < N Z T ) =* 0 7. 1 L F. ( N T A B T , 1 >- 0 Z 1 LE< 1 » 1 ) 

C ZFRAOS UP TO THROAT- 

POWER 1 = 1 5 
FNZ—FLOAT (NZT > 

DO 500 IZ=l,MZT-l 

ZFRAC (MZT-1 Z > = ' l . - ( FLOAT ( I Z ) /FNZ ) **PC*WER i. ) 4 IFF LAC (NZT \ 

500 CO NT INUE 

C ZFRACS DM3, OF THROAT — 


C 


FZN07— FLOAT ( M7--MZT ) 

Z NO Z« ZFRAC ( NZ ? -ZFRAC ( NZT > 

00 501 .17 “NZT -i-l , MZ - 1 

501 ZFRAC ( 1 1 > = ZFRAC ( NZ T ) + ( FLOAT < IZ-MZT ) /F'ZNOZ > **P0MEF2«ZN0Z 
-ADD AN ADDITIONAL CELL DNS. OF EXIT PLANE. 

ZCLAST-2. 4 ( ZFRAC < NZ > -ZFRAC ( NZ-1 > ) 

ZFRAC ( NZ + 1 ) - ZFRAC ( NZ ) -t-Z CLAST . 

NZ-NZ + 1. 

WRITE (6, 6723) (ZFRAC(IZ), IZ = 1 ,NZ>. 

723 FORMAT ( " ZFRAC" / ( 1 F3E .1 i . 3 > ) 


C GROUP 6 . MOVING GRID OR DISTORTED (BODY-FITTED) GRID s 

C - DISTORTED (BODY-FITTED) GRID : 

C MD, 

C XE ( 1 -MYF 1 , 1 — N 7. P 1 , l-MDX ( NYP1*NZP1*MD) *1 . 0>, 

C X W ( 1 -NY PI, 1 -N Z P .1 , 1 -MD ) . 

C VN ( 1. NX PI , 1-M7P1 , l-MDX (MXP1 #NZF1 *MD > #1 . 0>, 

C YS ( 1 -NXP 1 - 1 -NIP 1 , 1 -ND ) , 

C Z H ( 1 -MX PI , 1-MYP1 , l-MD) < ( NXPli'-NYF 1 *ND)*1 . 0>, 

C ZL ( 1 NXP 1 , 1 —NY' PI , 1 -MD) , STOPS' A ( 1-6X6*. F. > , STORWD ( 1 -6 ) <6*.-F . > , 

C STORP< . F . > , STORPEC , F . > , STOPPNC . F . > , STORPH< . F . > , F R T E4"-'C< , F . > , 

C DARCY 
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NOTES FOR VELOC I TV -F I ELD PR I HTO'JT 

( 1 ) THE STANDARD VELOCITY-FIELD PRINTOUT IS FOR THE 
VELOCITY RESOLU1 ES AND IS ACTIVATED IN THE USUAL 
WAY . 

(2) AM ADDITIONAL OPTION EXISTS FOR PR I NT 7 NO THE 
CARTESIAN VELOC I TV -COMPONENTS WHICH NAY BE 
ACTIVATED BY SETTING THE FOLLOWING LOO I CALS ' 

A) S10VARNJ2)-. T. FOR U-COMPONEN T 

B ) STOVAR ( V2 > - . I . FOR V — C 0 M FOMENT' 

C) S TOVAR (W2>=. T. FOR W— COMPONENT 
ftNAF’N 1. NO ! ! ! ! ! ! 

WHEN USING DECS PLEASE NOTE THAT STOVAR ( H3 > , STOVAR ( C4 ) 
AND STOVAR (21 ) ARE UNAVAILABLE. 
tHISER SETS ND (NUMBER OF SUB-DOMAINS) HERE : 


BFC= - T RUE . 
NU-- 1 


TOR 

:sa ( 

3 ) - . 

TRUE. . 

TOR 

WD ( 

3 ) - . 

TRUE. 

TOR 

.3 A ( 

4 ) = . 

TRUE , 

TOR 

SA ( 

5 ) . 

TRUE 

1 OP. 


6 ) = . 

TRUE, 

TuC 

RTF- 

.FAL 

SE . 


CALI.. DOHA I N ( 1,1, NX , 1 , NY' , 1 , MZ ) 

1 7. TAB“ 1 

7. ~ 07 ZLE ( I ZTAB- 1 ) 

Y—OZZLE ( I ZTAB, 2) 

DO SO). J Z •-=!, NIP 1 
IF(IZ.EQ.l) THEN 
YNN=Y 
GO JO 600 

ELSE IF(IZ.EQ.NZ) THEN ■ 

D I WwQZZL E ( NT AB , 2 ) -YNN ■ 

YNN--OZZLE ! M TAB , 2) 

GO TO 600 

ELSE IF ( T7.FQ.NZP1 ) THEN 
YNN~ YNN +D I VV 
GO TO 600 
END IF 

7 Z 7. * Z FRAC < I 7 - 1 ) - Z FRAC ( N Z T ) 

.104 I ZTAB-* I ZTAB+ 1 

I F ( IZTAB. 0 T . N I A B ) 0 0 T 0 6 0 1 
ZLAST — Z 
' YLAST«Y 

Z— OZ ZLE ( IZTAB, 1 > 

Y-OZ ZLE ( IZTAB, 2) 

F. F ( Z + . 00000 1 . LT . Z Z Z ) 00 TO 1 04 
YNN==YLAST+ ( Z ZZ -ZLAST ) * ( Y-YLAST ) / ( Z-ZLAS T > 

600 DO 602 I. X = 1 , NXF- 1 

YN ( I X , I Z , 1 ) ='r NN-*-r THROAT 
602 VS ( I X , I Z , 1 ) -0 . 

601 CONTINUE 

DO 605 I Y-- 1 , NVF 1 
DO 604 I Z — 1 , N Z P 1 
I F ( T.Y.EQ. 1. ) THEN 
X XE=0 . 

ELSE 

X X E-- .. 0 1 * YFR AC ( IY-1 ) *YN ( I , IZ, 1 ) 

END IF 

X E ( I Y , I Z , 1 ) -• X X E 

604 Xl-J( IV, IZ, 1 ) - - X X E 
DO 606 I X= 1 , NXF 1. 

ZH (IX, I V , 1 ) = ZFRAC ( NZ ) # THRO AT 
606 ZL ( IX, IV, 1 ) — 0, 

605 CONT I Ml JE 
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0 

f 
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r 
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I 
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i' 

i 
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i 
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i 

i 


i* 

t 


GROUP 

S: . I? 

EP 

EHP 

EM 

T 

i • 

AR I ABLE 

\Z* 

TO 

BE SQL 

VED 

rr* 

R OR 

STORED 

SOLVAE 

( 1. - 

,**" j 


■=; j'- 

. F 


> , ST OVA 

F. < 

1 - 2 

5 > < 23 -^ 

., F 

, r 

t jNO: :). ( 

1 ) ,;' 4 :> 

USE* FO 

LLO 

WI 

MG* 

NA 

ME 

n 

T NTEGE 

R 3 

FO 

R AP R A 

•/ F 

LEf'i 

EM TS 

1 - 20 : 

P 1. . FP , 

U :!. , 

U2 

? V 1 

\ i 

2 ■» 

H 

i. . us , r-i I 

7 M 

2 7 F 

S 7 KE 7 E 

f : ' , F! 

1 , F 

2 ,. H 3 •, 

C 1. , C2 • 02 

'*:*( JL.'v AR 

(FI 

) -- 

. ]'R 

UE 
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SDL VAR 

( PP 

\ “ 

. r F 

IJE 

, 










SOL VAR 

( V 1. 

) 

. TR 

UE 











SOLVAE 

( H 1 

) = 

. TP 

UE 











SOL VAR 

( C 1 

) •*- 

. TR 

UE 

„ 










SUL.. VAR 

( KE 

) 

. TP 

UE 

T 










SOL VAR 

< EP 

) 

. TR 

UE 











SOI... VAR 

(HI 

) ~ 

. TF 

IJE 

, 










S TOVAR 

( 23 

) = 

. TR 

IJE 

a 










S 1 OVAR 

< H2 

) = 

. TF 

UE 

, 










S r 07 A E 

( Uf 


. TR 

UE 

» 










stovae 

{ v? 

) ::rt 

. TP 

UE 











s rovAR 

( M 2 

) ^ 

. IP 

UE 











EX ITT'- 

ST 

nr; 

AGE 

F 

0F ; 


C OMTFtAV 

AK 

T AO 

T VEL .0 

C I T 

V 0. 


ENTS 

S TOVAR 

( ,t 6 


. TR 

HE 











■::« f UVAR 

( 2 1 

) = 

. TF 

'IJE 

• 










GROUP 

O 

VA 

R I A 

i Bl- 

E 

!.. 

ABELS : 








TITLE ( 

l --2 

5 ) 

<2F 

ip 1 

. 2 

H 

FP, 2 HU 1 


HU 2 

: . 2H v i , 

2 HV 

-•i 7 id 

'HW 1 7 2 

FIN 2 ? 2 HE J. 

2 HR 2 , 2 

HRS 

? 2 

HKE 

7 .,-2 

HE 

P 

, 2 HH 1 7 2 

HH 

*“.! 

HH 3 , 2 H 

C 1 7 

2 HC 

■“* 


2 HC 3 , 2 

HC 4 

■» 2 

! IPX 

v 2 

HR 

:V 

» 2 HRZ ’ 

2 ; '- 

4 H« 

• 





T I TLE ( 

2 ) 

= 4 

HRH 

0 1 











r I TLE ( 

H2 ) 

M 

HTE 

.rip 











T I TLE ( 

23 ) 

-4 

HCE 

ER 











GROUP’ 

1 0 

PR 

OPE 

R ! 

IE 

2 









I RHO 1 < 

: i > , 

IP 

HO 2 

:<i 

7 

R 

1:10 1 < :!. , 0 

’> , 

RHC 

12 < 1 . 0> 





AR HO 1 < 

1 . 0 

> 7 

BRH 

!0 1 

< 1 


0> , CRI 4 Q 

i <: 

J. . 0 

l 





I EIIU 3 . < 

1 > , 

EM 

U 1 < 

: i .. 

0 “v 


EI'iULAM< 

i . 

E -1 

0 > 





I H® AT . 

H 1 S 

AT 

, H2 

3 A 

r. 

P 

3 A TEX 2 1 

. 0 







S ). GMA ( 

1 -2 

3 ) 

< 1 . 

( ) 

£ , 

o 

7 1. • * 1. , E 

1 0 

7 1. . 

7 1 . E 1. 0 

7 .1. . 

7 1 . 

E l. On 


4 * 1 . 0 , 

1 . 3 

14 

, .1 . 

0 7 

i . 

F 

.1.0? 10*1 

. 0 








SIGMA (24)-. 7 
SIGMA ( 1 4 ) - . 88 
EMULAM-1 * 02E-4 

EMU.1 “EMULAM . ■ 

T RHO .1. =“ J. 

MOLAR MUMPERS ; SC<1> FOR H20 SC < 2 > FOR H2 SC < S > FOR 

MASS FRACTIONS: SMO ( 1 ) FOR H20 SMO < 2 > FOR H2 SM0(3> FOR 

FMI X~6. 054S51 

SMO < i >”-0. ' 

SMO < 2 > - 1 , / (.1 . +FMIX > 

SMO ( 3 ) - 1 . -SMO ( 2 ) 

ENT HMX~— 55. 7768E C 5 

WR I T E ( 6 * 6783 ) SMO ( !. > , SMO < 2 > . SMO < 3 > 

783 FORMAT ( M H20 ( 0 ) H2 \ 0 ) 02 ( 0 > " , 1 P3E 1 2 . T > 


0 ;> 

i j 


E-3 



) 


i.„:- 


SELtC'r LUMPUS »; i 1 Of i f 0 Nh I'CH AljlnBAnC TEMPERA'! LJPE 

DO 1 777 I TRY- i - 1 

El 2 0 A 0 = , O 4 3 4 1 7 5 + 1 E - 7 # F L. 0 A T ( I T R V - 1. ) 

SMB ( 1 ) =8 MO ( 1 ) 4- H 20 AD 
SMB < 2 > «SM0 ( 2 > -H20AD/ ? , 

SMB ( 3 ) ™*SM0 ( 3 } --H20AD* 3 . /? . 

C MASS FRACTION /MOL NT OR MOLES/UMIT MASS, 

SC ( J ) —SMB ( 1 ) / 1 S - 


TEMPO™ 300 , 

CALL EM rHAL < 7 EMF'O , EMTH , CP DR , SC , 3 > 

H H H - E N T H # 1 E MR 0 # R G A S 
E E E - C P D R * R G A S * ‘ !“ E M P 0 

UR I TE ( 6 , 6749 > SC ( 1 ) ■* 1 S . ■ SC < 2 ) * 2 . SC ( 3 ) *32 
.! HHH , EEE , HHH-EEE 

6749 FORMAT ( 11 ASSUMED INLET COMPOSITION AND IE 

1 " M ( H20 ) = : 1 PE .1 2 . 4 , " M ( hi 2 > - " , .1 PE 1 2 * 4 , M 

2 " TO CP H E HO” ? 1P6E12.4) 

• MR J 7 E ( 6 , .1 5 1 9 ) H20 A D , EM F I IMX , HHh ! 

1. 5 1 9 FORMAT ( H H 20 AD E X EE ” , j. P3E 1 2 . 4 ) 

1777 COM I I HUE 


L , TEMPO ’ CPDR^-RGAS, 

.! 'i r’E PAT URE M / 

M(02)= " , IPE.12,4/ 


r; 


p lh . #SC ( 2 > 

C M ;i XTUPF AFTER; COMBOS I 1 ON 

SC ( 1 ) = SC ( 1 ) +2 , *3C < 3 :■ 

SC ( Si ) =o « 

•FMFB— 2 - *SC ( 2 > 

TTEMFR- TEMPO 

CALL TEMPER ‘ HHH , TTEMPR , TEMP - CFDP , ROAS - SC < SO 
1.0(2 I TE ( 6 , 6709 > T TEMPR ? TEMP - HHH 
6709 FORMAT ( "TEMPER FOUND TO T H" - 1 P3E 1. 2 . 4 > 

C A L L E N I I* \ A L ( T E T i P , E N TH , C P D R 7 3 C , - 3 ) 

HR 17 E ( 6 7 675 1. > TEMP , EMTH , EM TH*RGAS aTEM-* 

675 1. FORMAT ( " V E NT H M , i P3E ! 2 . 4 ) 

Hr- I f E ( 6 7 4 750 ) SC < ) * 1. 8 . , SC ( 2 ) *2 . ? SC ( 3 ) #32 . 


4750 FORMAT ( " FULLY BURNED COMPOSITION AND ADIABATIC 1 EM SELECT f ON' 
j. / M M ( H20B > =■ " • 1 PE 1. 2 . 4 , " M i M2F ) = * l ' - 1 PE 1. 2 . 4 , " M ( 02P ) - " , :! PE !. 2 . -1 

C - M E A N M 0 L E C L 1 !... A R !■ * E I G M 1 * . A • • i n V m m A | p. ■. 1 • r 

H TMOL-- ( SC < 1 '* * 1. S „ +SC ( 2 ) #2 .. ) ( c -» • 1 i-SC f 2 ) '• 

hi 23 AT ~ R 0 AS *-Ef- 1 1 1 1 

C CHAMBER AREA 

RCHAM= 2267 

A I NLET=RCHAM* *2 *3 . 1 4 1 59 
FLOW A 1 -AMASF/A I MLET 

C A/A* AT INLET 

A AT ~0Z ZLE ( 1 ? 2 ) **2 
AM =0 . 

■ CALL. MSOLV ( GA , 0 , A AT , AM , S , TTT , RTF , RTF; , VRT.» QRT > . 

T TOT ~TEMF : 

PTOT-FL.ONA 1 *SQRT ( TTOT / WTMOl.) /QRT 
P TOT-2935. 7 * 6994.757 
WFAC=SQR1 ( TTOT /NT MOL ) 

C INLET VELOCITY 


WIN=VRT*WFAC 

l-JR LIE ( 6 , 674 4 ) W I N , TTOT , PTOT 
6744 FORMAT ( " N I N = " , 1PE 1 2 .4," TTOT F TO T " , .1 P2E 1 2.4) 

C TOTAL ENTHALPY 

H1SAT-TT0 T*H2SAT 

WR I TE < 6 , 6497 ) A AT , AM , TTT , FTP , RTR , VRT , QRT' , TTOT , H I 
1 UTMOL, <SC( II ) , I J=1 ,3) 

6497 FORMAT ( "A AT AM TTT FTP 

1 1 P7E 1 2 . 4/ " TTOT HI SAT H2SAT 

2 1P7E12.4) 


SAT r H23AT., 


QRT" / 
SCI - 2-3" / 


RTR VRT 

. WMOL 



c OUTLET PRESSURE 

A A DO Z ZLE ( N TAB , 2 ) **2 
AM “2 . 

CALI.. fISOLV ( GA , 1 , A AT , AM - S , TTT , FT P , RTF: , VRT , QRT ) 

PQIJT-P TOT/FTP 

WE I TE ( 6 , •'•■297 ) A ATT AM , TTT , FTP - RTR , VRT , POUT 
6297 FORMAT ( "AA'T M TTT .RTF RTR VRT POUT" . 

.1. / 1. P7E 1. 2 . 4 > 


C GROUP 1 ! I NT ER-FHASE TRANSFER PROCESSES s 

C I CF [ P , CF I PS , I MOOT , CMDOT , CA 1 I < 1 . E6> , CA2 1 < 1 . E6> 

C 

c GROUP 12 SPECIAL SOURCES s 

C I SPCSO < .1 -25 ) , AGR A V X , AGRAVY , AGRA V Z , ABUOY , HREF 

C r GROUP 13 INITIAL FIELDS : 

C F 1 1 NIT ( 1 -25 ) <25# 1 . E - 1 C» 

FIIN.n(P.1 ) =10101 . 

F 1 1 M J. T < U 1 ) = 1. 0 1 0 1 . 

I F ( SOL VAR ( C 1 > ) F 1 1 N I T ( C 1 ) = .1. 0 i 0 1 . 

FT I MIT 0-11 ) =H 1 SAT 
FI III] TOT 2) =10.1 01 , 

C-- FOR INITIAL FIELDS FOR K AND E SEE BOTTOM OF CM. 1 * 

C ‘ GROUP 14 BOUNDARY/ INTERNAL CONDITIONS 3 

C I LOOP 1,11. CORN , X CYCLE< . F . > , P BAR , REG I ON ( 1-1 0 ) < 1 0# . T . > 

C #N.B. ALL 10 REGIONS ARE DEFAULTED ..TRUE,, THE USER SHOULD 
C SET REGION! I >~. FALSE. FOR UNUSED REGIONS I •" . 

DO 140 IP~ 1,10 
1 40 RE G I ON ! I R ) - . F AL SE . 

C GROUP 15 TO 243 REGIONS 1 TO 10 

C- ONLY THOSE REGIONS ARE ACTIVE WHICH ARE SPECIFIED BY THE 

C USER, PREFERABLY BY WAY' OF : - 

C CALL PLACE ( IREGM* TYPE, IXR, I XL, I YF , IVL, IZF , IZL) C 

C CAL L COVAL f I REG I -I , VARBLE , COEFF , VALUE ) 

CALL PLACE ( 1 , LOW, 1 - 1,1, 20, 1 , 1. ) 

C CAL L COVAL < 1 , M 1 , F I X FLU , FL OWA .1 ) 

CALL COVAL! 1, Ml, ,05 •• FTOT ) 

I F ( SOL VAR < C 1 > ) W I N«W I N#TEMPO/ 1" TO T 
CALL COVAL ( 1 , W .1. , ONLY! IS , W I M ) 

.1. F ( SOL VAR ( C 1 ) > CALL COVAL ( 1 , 0 1 , ONI.. VMS , FMUB ) 

I F ( SOLVAR (HD) CALL COVAL ( 1 , H 1 , QNLYMS , H 1 SAT ) 

CAL L PL ACE ( 2 , H I GH ,1,1,1, 20 , N Z , M Z ) 

CALL COVAL ( 2 , M 1 , . 05 , POUT -2000 . ) 

IF ( SOLVAR ( KE ) ) THEN 
TIJRBLV-.l 

SLEMTH=. 02#THR0AT#(0ZZLE( 1 ,2) > ' ’ * 

EK I N“ . 5# < W I M#TURPLV ) **2 
EP I N= . 1. 6 4 3 # E K I N # # 1 . 5/SLEMTH 
WR I TE < 6 , 543 1 ) EK I N , EP I M 
543 1 FGPMArC" EKIM ERIN ",1P2E12.4) 

C CALL PLACE ( 3 , NOR TH ,1,1, NY , NY , 1 , NZ ) 

C CALL COVAL ( 3 , N 1 , WALL , 0 . ) 

C H.1 VALUE WILL FE OVERRI TEN IN GROUND 

C CALL COVAL ( 3, HI , HALL, l.) 

CAL L COVAL ( 1 , KE , 0! IL VMS , EK I M ) 

CAL L COVAL ( 1 , EP , ONLY MS , EP I N ) 

C INITIAL FIELDS FOR K AMD E 

FI IM1T(KE>=EKIH 
F 1 1 N I T ( EP > =EF I M 



c GROUP 25 GROUND STATION s 

GROST A< . F . > - M AML ST < , F , > 

C SNA MLS T ACTIVATES NAMELIST IN GROUND. 

C GROUP 26 SOLUTION TYPE AND RELATED PARAMETERS : 

C WHOLEPC . F . > , SUSP'S T< F , > , OONAPC 2 , F . > 

LOGIC ( 87 ) SET TO . TRUE. FOR COMPRESSIBLE FLOPS 

LOGIC ( 87 > TRUE . 

C GROUP 27 SNEER AND ITERATION NUMBERS' s 

C FSNEEPC 1 > , LSNEEPC 1 >. LI THYLX 1 > , LI TC< 1 > , L I TKE C l >, LITHC1>, 

C L I V ER ( 1 -25 > <9 *1,-1, 15*1 > 

C I V E L F < I. > , N V E L < 1 > , I V E I... l_ < 1 0 0 0 0 > , 

C KEFC 1 > , NK.EC 1 > , I KELT. 1 0000> , 

C I ENTF< i > , HEMIC 1 > , I EMTLC J 0000> , 

C I Cr.!CF< .1 > , NCNCC 1 > i I CNCL< l 000 0> , 

C I F I 10 1. F< 1 > , MRHO 1 < 1 > , I RM0 1. L< 1 0000> - 

C I RH02FC .1 > , NRH02C 1 > ; I RH02L-C1 OOOO 

F SWEEP- 4 01 
LSNEEF=450 

r - g R o u r : ' 2 3 ■ r e r m t n a t i o n c r i t e r t a : 

C END I T ( 1 —25 ') <8 * ! . E- 1 0 , 0 . 5 - 1 5* 1 . E- 1. 0> 

C GROUP 2? RELAX ATI ON : 

C R L X F < 1 , > , R L X P X Y < t . > , R L. X P 7 < 1 , > , R L X R H 0 < 1 - > - R l . X M U T < 1 . • 

C DT PALS < 3—25 ) <23* 1 . E 1. 0> 

RLXP=. 6 


RLXPZ-.6 
RLXf-'X Y— - 6 

DTFAL.S ( o 1 ) . E3 8-THRO A F 

DTFAL.3 < V .1 ) =D FFALS ( W 1 ) 



GROUP 30 LIMITS : 

VR.HhXO, . ElO'T- V’ELM I f-K.-l . E10>, RHOMAX-C 1 F. 1. 0> , RMON I N< 1 
TKEMAXCl . E 1 0> » TKEM I M< 1 . E-10>, EMUNAXCl . E10>» EMUMH-K1 
E P ’SM"AX<1 . E 1 0> ’ EP S M I M< I . E - 1 0> • AMDTM X< 1 . E10>, AMDTHMC- 


E- 1 0. 

e-i o: 
. f ! o: 


GROUP 3 1 SLOW I NO DEV I CES : SLORHO C 1 . > , SLOENIK 1 
SLORHO— . 3 


C — GROUP. 32 PRINT-OUT OF VARIABLES : 

C PR I NT ( 1 -25 ) < . T . , , F . , 23* . T . > , SUBWGR< . F . > 

PR I NT (.23 ) = . TRUE . 

PRINT (2) =: TRUE. ' 

PR I NT ( H2 ) = . TRUE . 

I 7 MON -12 
I YMOM-“ 1 9 


GROUP 

-Ji 

Mi 

ON I TC 

JR P! 

R I NT -OUT 



- 

I XMON< 

1 > 

, I 

YMOM< 

■- 1 7 

I Z NONE 10 

•, NPRMOM 

I<1>, 

NPRMNTC 1 1 > 

GROUP 

34 

F 

I ELD 

PR I 

NT -OUT C 

OMTROL 

5 


NF'R I NT 

<; i 

oo : 

> , n it 

: ’R I N 

< 1 00 > J NX 

PR I N <1 > 

, HYP 

R I N< 1 > ; N Z PR I N< 1 0 - r 

I ZPRF< 

1 > 

, i 

STPRF < 1 > 

, IZPRLC 1 

OOOO t I 

STPR 

;L< 1 0000 > 

NUMCLS 

< 1 

o > 

, K 0 U 1 

'FT 





NPR I NT 

=2 

00 







NYFRIM 

--2 








KOUTPT 

— -- 

1 







GROUP 

35 

Ti 

-OBL. E 

CON 

TROL : 




TABLES 

. 

F. 

>, NT ABLE 

, MTABVR , 

L INTAB' 

NFRT 

AB,NMON, 

I TAB ( 1 

O 

) - 1 

■IT APT 

/R ( 1 

- 3 ) 
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7 


C CROUP' 36-38 ARE NOT DOCUMENTED IN THE INSTRUCTION 

C MANUAL AND ARE INTENDED FOR MAINTENANCE PURPOSES ONL «' 

C OROUI : ' 36 DEBUG PRINT-OUT SLAB AND TIME-STEP s 

C I ZPRK1>» I Z PP2< 1 > j I STPR l < J > , ISTPR2<1> 

C GROUP 37 DEBUG SHEEP AND SUBROUTINES s 

C MENU • KM A I N , K INDEX , KOEON, K INPUT , KSODAT - KCOMRF, KSORCE- 

C KSOI 1. , KS0LV2 - KSC." .VS , KCOMPP , KADJST , KRJ.J X , KSH I FT , KD T. t:: , 

C KCOMPUi KC0I1PV, KCOMPW, KCOMPR » KW ALL , K!JBRHO<— 1 > , KDBEXP. KDPHDT 

C KDBOEM 


OR 

OUR 


MON I TOR , 

TE 

ST 

, AND 

FL 

AG : 








HON I TR 

< , F 

. >,F!..AG< 

. F 


, TEST 

T n 

T , > , K 

FL 

.60 

<: 3. > 





END OF 

t J f\ 

|l H 

I NTENANC 

E~ 

ON 

LY SE 

CT 

ION 








GR 

UUP 

39 

ERROR AN 

D 

RE 

S .1 DUAL 

PR I NT 

— c 

H JT 

: 





I EERPC 

1 000> , RESRE 

.F ( 

3. , 

3-24 ) 


5- 3. . > 

,F 

; '.Ec 

MARC . F 

• S' ? 




RE 

S I D ( 

1 -2 

5 ) <2-* . F . 

T jC 

•3; y.. 

. T . > , 

KGUTPT 








RE 

SEE F 

• i ) 

-WIN^AIN 

*LE 

T* 

. 0 1 










GR 

;;njp 

40 

SPEC I AL 

DA 

TA 

: LOG I 

C ( I . . 

.1 C 

)) 7 

I NTOR ( 

.1. . . 

1 0 

) 7 F 

E ( 2 3. , . 3 

in... 

3 PCI. 

> , N 

I SRC 1. > , 1 

IRS 

F< 

3. > , SF : 

DA 

[AC - F 

»' -• 

;• 7 L 

3RD A ( 3. 

) , I 

SP 

DA ( 

J. > , Rpr-D 

IJ:? 

E FI 

EST 

10 EL.EI- 

EM 

TS 

OF A 

FR 

H T 3 L 

OF 

1 1 C 

« INT 

OR 

AM 

Fi 2 

i ST 

TO 

SOT 

H OF ARRAY 

RE 

F 

OR IF 

AM 

3FERR 

JU- 

■IG 

SPEC I A 

!... D 

AT 

A F 

'ROM 

3 A 

TELL 

ITE 

TO OR Oil 

HD 


BUT I 

F 

REQU I 

RE 

•NF. 

NT'S EX 

CEE 

D 

TH I 

'§ 

PR 

0 V I s 

I ON 

SET 3 FT 

lAT 

A 

- . T , 


AND L 

4. t 

1EH 

SION A 

RRA 

r !-;i 

LS 

PDA , 

IS 

PDA, 

RS 

PDA A.BOv 

F 

AN 

D IN 

OR 

OUND 

f\ f*" 
!- 1 

NEEDED, 

AND 

SET 

HERE 

ML 

SP * 

N I S 

F, NR3P 

TO 

D 

I MENS 

T 0 

N VAN. 

UE 

:S., 






GR 

OUP 

42 

RESTARTS 

; A 

ND 

DUMP 

C; 

: 8AV 

Er 

1c . 

F. >, RE 

31 E 

T C 

1, F . 

- k t r !i‘ ij 

SA 

■vEM-" 

, TR 

UE . 













I F 

( F3M 

ERR 

J3T.1 ) F 

:FS 

TR 

T= . TR 

UE 

m 









C X X X X X X 

cxxxxxx 


X X X X X X X X X X X X X X X X X 
X X X X X X X X X X X X X X X X X 


XXXXXXXXXXXXXXX USER SECTION 3 END 
v X x X X X X X X X X X X X X STANDARD SECT I ON 4 


ART 


C WRITE GENERAL DATA ON TO THE GUSJEt.DTA TAPE, ETC... 

I F ( SPDATA ) CALL WRTSPC ( LSPDA , ML3P , I SR DA , N I SP , RSPDA , MPSF 
1 1- ( BLOCK ) CALL WRTPOR ( RE , PN , PH * PC , NX . NY , NZ , I PLANE ) 
IF(BFC) CALL NR T DEC ( 14, NBFC, XE, XTL YM, VS, ZH, ZL, 

&ND, NX+1 , MY+1 , N Z +■ 1 , NZ , PR TBFC ) 

C OLD PRACTICES RETAINED FOR REFERENCE J 

C IF (SPDATA) CALL SPCDATURIJN) 

C IF (BLOCK) CALL FORD AT ( I RUN > ... 

IF (GRAPHS) CALL SORT (I RUN) 

IF(REST'RT) GO TO 902 
DO 901 I ND V AR= 1 , 25 

I F ( I F I X ( F )' I N IT ( I NDvAR ) +0 . I ) , ME . 1 0 1 0 .1 ) 00 TO 90 3. 

NR I T E < 6 , 6742 ) FMUB , FMFB , H .1. SAT 
6742 FORMAT (" FMUB BEFORE CALL FLDDAT ", 1 P3E 12.4) 

CALL FLDDAT ( I RUN) 

GO TO 902 
90.1 CONTINUE 
902 CALL DATA I 0 ( NET , 1 0 ) 

I F < MON I TR ) CALL DATA 1 0 ( WET , -6 ) 

999 CONTINUE 
STOP 
END 



i 
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SUBFOUT J. ME T EMF'F.R ( HSTAT - TO , T , CP DR , RGAS , SC , M3C > 

C SUB ITERATIVE CALCULATION OF 1 EMP'ERA’I URE 

DIMENSION SC (NSC > 

DATA N I TER , DTO 7 TM I N/ 1 2 , 50 . , 1 00 . / 

DT-DTO 
TENR“ T ( ' 

CALL EM THAL ( TEMP , FHH , CPDP , SC 7 NSC > 

ENT hi — H H ! i *■ R 0 A S T FIT- . 

I F i HST AT , LT . EN r H ) DT=-DT 
TEMPI..- I'EMP 

C NR T TE ( 6 7 6505 ) TO 7 EMTH 7 MS TAT , P0A3 - SC ( 1 ) 7 SC ( 2 > 7 SC ( 3 ) 

C6505 FORMAT < " TO E H3 RG SCMP7E12.4) 

TEMP -TEMP+DT 
I TER-0 

1.0 ENTHL-EM TH 
ITER=ITER+ :l 

CALL EN THAI... ( TEMP 7 HHH , CROP , SC 7 NSC ) 

EN TH «HHH*RGAS*TEMP 

RENT H~ < 1 1ST AT -ENT ML ) / ( ( E NTH -ENT HI. > ' - 1 - 1 . E-9 ) 

C HP I TE ( 6 7 6500 ) I T ER , T EMP , ENT H 7 EMTHL 7 HST AT , REMTH 

06500 FORMAT ( " FT T E EL H3 RE" 7 12- IP 5E 12. 4 ) 

I F ( APS ( E NTH -EMTHL > .LT. . 00 .1 *A8S < EMTH > ) REMTH- .! . 

TEMP 1 -TEMPL+ ( TEMP--T EM PL ■’ sRENTH 
TEMPI = AM AT I ( TEMPI 7 . 5 -K-TEHP 7 TM I N ' 

T EH P 1. -AM I M 1. ( TEMP' 1 7 1. . 5» TEMP 7 5000 . > 

TEMPI... - Tf-TTP 
TEMP-T EMP 1. 

AR-ABS ( REMTH > 

IF ( (AR.GT. 1.005 .OR. AR.LT. .995) . AMD. ITER. LT. NITER) GO TQ 10 

l‘=TEMF 

RETURN 

El ID 


SUBROIJ T 1 1 IE EN F HAL. ( TEMP’ , 
D 1 MENS r ON SC ( NSC ) , 7 S ( 7 , 
D AT A 7. S / 2 . 7 1. 7 , 2 . 9 45E- 3 , 
3. 7 6 . 63 1 , 4 . 07 O0E-3 , 

2 , -3 . 028E4 , -3 . 227E- 1 
3 - 7 3 .1,5.1 1 2E-4 , 5 , 264 E-S 
4 7 - 1 , 963 - 3 . 057 7 2 . 667E— 3 


hISUM 7 CP SUM 7 SC 7 NSC > 

; -I > 

-- S . 0 2 2 E - 7 , 1 . 0 2 3 E - 1 0 , -4. 8 4 7 E - 1 5 , 
4.1 5 2 1 — 6 7 — 2 . 9 >6 4 b. — V 7 S . O 7 b. — 1 3 

7 -3 . 49 1 E- 1 1 , 3 . 695E - 1 5 , -8 . 77--P2 
7 -5 . 8 1 E-6 , 5 . 52 1 E-9 , - 1 . 8 1 2E- 1 2 




6 7 3 . 622 7 7 .. 362E- A . — 1 . 965 E— 7 7 3 . 620EI— 1 1 7 —2 . 

7 7 .615? 3 . 626 7 • - 1 . 87SE— 3 7 7 . 055 E— 6- 7 —6 . 7 64 E 

8 , -1. . 048E3, 4. 305/ 

NS -NSC 

I F ( SC ( 3 ) . LT . . 00000 1 ) N8=NSC- 1 
K-- 1 

I F C TEMP . LT , 1 000 . 1 K=-*2 
TEN P 2 -• T E M P «• T E MR 
HSUM=0. 

CFSUM=0 . 


5E-.1 5,-1. . 202E3 
, 2. 156E-12 


on 6 I S= 1 , MS 
C F .1 = Z S ( 1 , K, IS) 

CP2-ZS ( 2 , K , I S ) #TEMP 

CF3=ZS ( 3 , K , I 3 ) * TEMP 2 

.CP4-ZS ( 4, K, IS ) *TEMF’2#TEHP 

CP5-ZS ( 5 , K , I S ) <- rEMP2*TEMP2 

CFSUN-CPSUM+SC ( I S > * ( CP 1 +CP2+CP3+CP4 + CP5 ) 

6 HSUM --HSUM+- 

1 S C ( IS ) » ( C P 1 + . 5 # C P 2 + . 3 3 3 3 3 * C F 3 + . 2 5 * C P 4 + . 2 C P 5 + Z S ( 6 , K, [ S ) 
RETURN 
END 


9 9 1 P.4 


•TEMP ) 
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TTT, p'[ P, RTR, VR OR F > 


SUBROUT I ME MSOL V ( GA - KS , AAT , A? 1 , S , 
DA T A R / c : -.:‘OG » £•/ 

G1“GA™ 1. 

02=0 .i -- - 5 
G3=- 1 , /G 1 
G*4 : -Gh/ U 1 
G5--GA+ 1 
06-05/ (2. *Q1 ) 

NITER=0 

1 F ( KS , G r • ) U U T u 

C SUBSONIC — 

T F ( AM , Or . , 9999 > AM«0 . 

1 AML -AN 

AM- < ( 2 - +0 1 *AM*AM ) /05 > **06/ AAT 

M I - TER-MITE'R+l 

IF (UTTER, OT, 100) GO TO 99 

I F ( APS < AM -AML ) / AN . L T „ . 00 ,1 ) THEN 

S — AM-&AM 

GO TO .10 

END IF 

go ro i 

e SUPERSONIC — 

9 G 7 = 1 / i.i 6 

I F ( AM . LT , I . 000 1 ) AM- 1 . 000 1 

2 S =- < ( A A T * A M > # * G 7 - 2 . / G 5 } * 0 5 / G .1 
AML —AM 

AM-SQRT(S) 

NITER-MITER+ 1 
I F ( N J TER O f . .1 00 ) GO TO 99 
I F < ABS ( AM— AML ) /AM . LT . . 00 1 ) GO TO 
GO TO 2 

99 WE I TE < 6 , 999 ) AM . AML , AAT , KS 
999 FORMAT! -* TOO MANY ITERATIONS 
10 TTT - 1 , +02*3 
PIP— TTT 4*04 
ETR-TTT ft #03 
' VRT-SQRT ( GA*R/7 TT ) -w-AM 
UR I — yt*!R 1 ( UA/ R ) 4 AM / 1 T ) 4 -4 U 0 
RETURN 
END 


1 0 

,3E10„4 f 14) 
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SUBROU f I MR FLOP AT ( I RUM > 


C FIELDS IS USED TO SPECIFY MOM UNIFORM (H.TTIi 

C EARTH SETS UNIFORM INITIAL FIELDS 10 l-UMIT 

C WHEN F I M I I T ( MPH I ? = 1 0 1 0 .1 . 0 W1 \ I Cl I [ 8 THE 8 .1 ON* 

C READS THE MOM-UNIFORM FIELDS TAPE THE COMTEt 

C ARE SET HERE. 

C IT IS ESSENTIAL TO PROVIDE SETTINGS FOR ALL 

C FOR WHICH FI I NIT (MPH I > HAS BEEN SET TO 10101 

C ORDER fl U 9 T FOLLOW THE STANDARD ORDER Of 

C NOTE: EARTH PRINT'S OUT THE INITIAL FIELDS II 

r: 

€:■**■* PLEASE MOTE THAT THE PHI ARRAY l-i 1.1 S T BE 

C PH I < .[ V , I X >==..„ , NOT AS PH J (IX, J Y ) = . . . 

C ALSO, YOU M U S T DIMENSION PHI TO THE MA5 

C DIMENSIONS USED IN THE RUNS DEFINED IN SAIL !• 

C DIMENSION PHI (NY MAX., NX MAX.) 


CHAPTER I PREL I M I WAR I ES 

1 1 MCI. U DE 9 , C NNGUSS I . F T H / 0 
T» I MCLUDE 9 , GU3SEQU I , FTN/O 
TIN C L U D E 9 , N 0 Z B F 1 . F T N 

C ’ USER -EQUIVALENCES' 

EQUIVALENCE < FMFB » RE ( 4 ) ) 

D I MENS I ON PH I ( 20 , 1 ) , FDTAFE ( 3 ) 

INTEGER FT APE 
LOGICAL FIRST 

I NT EGER F 1 , PP , U 1 , U2 , V 1 , V2 , W 1 , W2 , R 1 , R2 , R3 , EP , 
5iC3 , 04 

DATA F 1 , PP , U 1 , U2 , V 1 , V2 , W 1 , N2 , R 1 , R 2 , PS , KE , EP - 
SC: 3 , C4 / 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 1 0 ,11,1 2 , 1 3 , 1 4 , 1 5 , 1 <■ 
DATA F I PST , FT APE , FDTAPE/ . TRUE . , 9 , 4HF I EL , 4HDS 
I F ( F I PST > CALL TAPES ( FT APE , FDTAPE , 3 , 1 , 4* 
FIRSTS. FALSE. 

THRTS=THR0AT**2 


CHAPTER 


USER CHAPTER FOR 


SETTING INITIAL FIELDS 


NNZ=NZ 


IF ( PARAB ) NNZ-1 
DO 2001 iz=i;nmz 
DO 2000 MPH 1=1, 25 

I F ( MPH I . EQ . 2 ) GO TO 2000 
I F ( I F I .X ( F 1 1 N IT ( MPH I ) +0. 1 ) . ME . 1 0 1 0 J. ) GO 1 
I F < . NOT . STOVAR < MPH I > . AND . NOT . SOL VAR ( MF1 
SECT I ON STARTS HERE*#**##****#* 
C FOR EACH VARIABLE, MVAR, FOR WHICH F II NIT (MV 

C 10101.0 IN BLOCK DATA, THE USER MUST INSERT 

C STATEMENTS FOR SETTING HIS WANTED NON-UN I FOP 

C THE FORTRAN-STATEMENT MODEL TO BE FOLl.OWED E 

C EACH MVAR IS AS FOLLOWS... 


C 


4. FIELDS.. 

(MPH I 1 , EXCEPT 
•iL THAT EARTH 
•ITS OF WHICH 

THOSE MPH I ■'?. 

I . THE SETTING 
7 BLOCK DATA. 

7 K0UTPT=-1 .. 

SET AS 

■( IMUN 
[T, VIZ. 


: HI , H2 , HS, Cl , C2 • 

HI - H2 , HS, Cl . C 2 - 
, 17 , 18 , 13 , 20 / 
LD, 2HTA/ 
i- NX-S- NY ) 


0 2000 

II ) ) GO TO 2000 

'AR) IS SET TO 
HERE FORTRAN 
:M INITIAL FIELD. 
V THE USER FOR 


E - 10 



TTT , FTP , RTR - VET , GET ) 


IF ( Mr'H ] . WE . P 1 ) GO TO 2 
AA r-YN (1,11+1,1.) **2 /THRTS 

IF< I7,LT.M7T> KS-0 
CALL USOLV ( OA . KB , AAT , AW - S , 

PSThT=PT0T/F7P 

' r s' r f \ t t t o t / ' r t t 
I'JR I TE ( 6 , 6782 ) IZ, AAT , AM , TTT - RTF - PS TA F - T9 TAT , VRTTH-IF AC 
6782 FORMAT < "MSOLV CALL TM FLDDAT IZ = M2, " AAT All “ . .1 P2E J 2. -1 - •' 
1 " TT FT P T W ", 1P5E12.4) 

DO 201 I X = .1 , NX 
DO 20 1 I Y= 1 , MY 
201 PHI < IY, I X ) ~ PST AT 

HR I. TE < FTAFE ) < ( PH I ( I Y , I X ) , 1 Y- 1 , NY ) , I X = 1 , MX ) 

205 CONTINUE 

C 

C INSERT HERE ADDITIONAL FORTRAN SEQUENCES, SIMILAR TO THAT 

C GIVEN ABOVE, FOR EACH M VAR FOR WHICH FI I NIT ( MVAR > - .1 0.1. 0 1. . 

I F ( MPH I . ME . N :). ) GO TO 305 
WVEL-VRT*WFAC 
DO 301 IX~1 , MX 
DO 301 f Y~ 1 . NY 
301 PHI ( I V, I X )= NVR! . 

MR I TE ( F TAPE ) ( ( PH I < I Y , 1 X ) • I Y-= 1 , NY ) , I X= l , NX > 

305 CONTINUE 

I F ( MPH I . ME . H2 ) GO TO 405 
DO 401 IX =1, NX 
DO 401 J.Y = l,MY 
401 PH IX IY, IX>=» T3TAT 

WR I TE ( FT APE ) ( (PHI ( IY, IX) , I Y 1. , MY ) , l X = .1 , NX ) 

405 CONTINUE 

I F ( MPH I , ME • C 1 ) GO TO 505 
COMP=FMFB 

IF < IZ. EG. 1 ) COMP--FMUB 
DO 501 IX --1,NX 
DO 501. I Y-~ I. , N Y 
501 PH I ( I Y , I X ) "COUP 

WRITE (FT APE) ( (PHI ( IY, IX), I Y-- 1 , NY ) , IX»1,MX> 

505 CONTINUE 

L" w # -s- -K- >!• iv w ttSTtLISER SECT I ON ENDS HERE -v - : /r it -z'- i> -r L 

2000 CONTINUE 

WR J. TE ( 6 , 679 .1 ) 1 1 , T ST AT , PST. AT , WVEL. , COMP 
679 1 FORMAT ( " I Z T P l-J MH2 " ,13, 1 F4E 1 2 - 4 ) 

2001 CONTINUE 
RETI IRN 
END 


C$D I RECT I VE **CMM8F t** 

COMMON/FOB /FI (5000) . - 

COMMON / C 1 8 / ND >' C I C / KOORD 

COMMON / C I D/ KDBGO , KOBOMF , KDBGCD, KDR1MD, KDBMFX - KDBCDT , KDBFCS 
COMMON /C I E/KDBOS , KDB 1 NS 
COMMON/ C I F / I GEM 

COMMON / CR A / X W ( 2 1 , 42 , 1 ) /CRB > XE ( 2 1 , 42, 1 ) 

F • / C R C / Y S ( 2 , 4 2 , 1 ) / C R D / V N C 2 , 4 2 , .1 ) 

& / C R E / Z L ( 2 , 2 1 , 1 ) / C R F / Z H ( 2 ,21,1 > 

■v / ij R i , / R i ij j\| / 1 . R H / D A R i ... 't 

COMMON / CLA /STORSA ( 6 ) , STORWD ( 6 ) , S rORP . STOPPE , 9 'f ORPM, 

5 3 TORPH , 3 TOR 1 , ST0R2 , ST0R3 , PRTCOM , PRTBF C , 3 TOD PN 

COMMON /CLC /'BFPLOT 

LOG I CAL STORP , 8T0RPE , STORPN , STORPH , STOP 1 , SI 0E2 , STOPS , 

?< STORSA , STORWD , PRTCOM , PRTPFC , BFPLOT - S VOCRN 

COMMON/NOZ/WFAC, PTOT, TTOT, GA, NZT, THROAT , FMUB 
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SUBF.nijr imp GROUND ( IRM, I CHAP, I8TP, ISMP- I ZED, IMDVAR) 

* T NCLUDE 9 , CMNOUSS I . F-'TM/G 

* I NCL UL'E 9 ? OUSSEQIJ I . F T'M/G 
C: * I NCUJOE 9 - NHL T ST . FTN / 0 

CX X X X X X X X X X X X X X X X X X X X X X x X X X X X X X X X X X X X X X X STANDARD SECT I ON ). S TAP I S 
C VARIABLE NAMES FDR USE JN GROUND: 

COMMON / T YRE / CELL , EAST , WEST, NORTH- SOUTH, HIGH, LON, VOLUME- NALL 
COMMON/ VAR/P 1 , PR, U.f - »J 2, Vj , V2- WI • M2- R] , R2, R3- 
&KE , EF' -Hi , H2, H3, Cl , C2, C3- C4- RX , RY- RZ , SI , 32 
COMMON /VAROLD/F 1. 0 , PPO , !...! 1 0 , IJ20 , V .1 0 , V20 - W 1 0 , N20 - R 1 0 , R20 , RSO r 
S<KEO, EPO, HI 0, H20, H30, CIO, C20- C30- C40, RXO, RYO, RZO , S10, S20 
COMMON/ VARLOW/P 1 L , PPL , U 1 L , U2L , V 3 L , V2L , W 1 L , W2L , R 1 L , R2L , RSL , 
&KEL , EPL , H I L , H2L , H3L , C 1. L , C2L , CSL , C4I.. , RXL , RYL , R7.L - 3 1. L , S2L 
COMMON / VARH I / P .1 H • PPH , U 1 H , I.J2H , V .1 H , V2H - IT .1 H , W2H - R .1 H , R2H , RSH , 
?<KEH , EPH , H 1. H , H2H , H3H.C1H- C2H - C3H , C4H , RX.H , RYH , RZH , S 1 H , S2! I 
COMMON / OMTRY / VOL , VOL U , AEhS T - AMuFl I H , AH I L-H , ALL 1 X , hMD v , AHD .i 
COMMON/ PROP/D 1 , D2 , D 1 DP , D2DP , HU 1 , HU ). LAM , EX CO , CFP , MD X , HST 1. , MS L 
COMMON /FRPOL.D/D J. 0 , D20 
COMMON / FRF LON / D 1. L , D2L , EX COL 
COMMON /'FRPH I /D3 H, I '2 IT , MtJ .1 H, EXCOH 
COMMON / V ARNX / XO , X»J- DXU- D X G 
COMMON 'VARMY/YG, YV- DYV, DYG, R, RV 
COMMON / VARNZ / ZG , ZW1 - DZW, DZO- W0R10 
COMMON /ODMSC I / X PLANE , YPLAME , ZPL AME ■ l T NO 
■ COMMON/ODHSCL/LSLAB, MSLAB, H3LAB, LAilMU 
REAL NORTH, LOW 

I NT EGER P 1 , PP , U 1 , U 2 , V 1 , V2 , W .1 , N2 , R 3. , R2 - RS - 
&EP , H ,1 , H2 , H3 , C .1. , 02 , C3 - C4 , R X , RY -RZ, S 3 , S? 

I NT EGER P 3. 0 , PPO , U 1. 0 , IJ20 , V 1. 0 , V2Q - W 1. 0 , W20 , R 1 0 > R20 - RSO - 

?<EFO, HI. 0, H20- H30, Cl 0, C20, C30 , C40- RXO, RYO, RZO, S.1 0 . 32 « ' 

I NT EGER P 1 L , PPL -01 L , U2L - V 1. L , V2L , W I. L , U2.L - R 1 L , R2L - RSL - 
&EPL , Hi. L, H2L,H3L , Cl L • C2L , CSL , C4L - RXL-RVL-RZL - 311. .• S2L 
I NT EGER P 1 H , PPH , U 1. H , IJ2H , V 1. H , V2H , W 3. H , II 2 H , R !. H . RSH - RSH - 
S<EPH , H 1 H , H2H , H3H , C 1 H , C2H , OSH , C4H - RXH , RYH , RZH , S 1 H - 32H 
I NT EGER VOL , VOLO , A EAST , ANORTH - AH 1 GH AEDX , AND'/ - AMOZ 
I NTEOER D 3. , D 1 DP , D2 , D2DP , EX CO , CFP , HST 1 , HST 2 
INTEGER DIO, D20- D 3 L , D2L , EXCOL , D 1.H- D2H, EXCOH 
I NTEOER XG , XU , DXU , DXG , YG , YV , DYV , DYG , R , RV - ZG , Z W I , DZW , 

S--DZG, WGRID 

I NTEOER X PLANE , YPLAME , Z FLAME 
L OG I CAL LSLA3 , MSLAB , HSLAB - LAMMU - LSPDA 
EQU I VALENCE < Ml , R 1 ) , ( M2 , R2 ) 

C SAIL IT -EQU I VALENT I RUN: 

EQU I VALENCE ( I RUN - I NT OR ( 1. 1. ) ) 

CX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X STANDARD SECT 1 ON 1 ENDS . 
CX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X USER SECT I ON 3. STAR TS : 

C ARRAYS ( DIMENSIONED NY, NX ) FOR USE WITH -ADD": 

D I MENS I ON CVAR ( 20 , 1 ) , WAR ( 20 - 3. ) , CM ( 20 , 1 ) , VM < 20 - 1 ) , Z EPO ( 20 , 1 » 

• C SPECIAL-DATA ARRAYS DIMENSIONED S: DIMENSION VALUES SET HERE; 

DIMENSION LSPDA ( 1 > , ISWA( 1 > , PSFDAI 1 ) 

C USER PLACES HIS VARIABLES, ARRAYS, EQUIVALENCES ETC. HERE. 

D I MENS I ON GV ( 20 , 1 ) , GW < 20 , 1. ) , OF ( 20 , 1 ) , GP I Z 1 < 20 , 1) , GRH ( 20 - t > - 

1 GTEMP ( 20 , 3 > j CFLH20, .1 ) , GZMODE ( 4 1 ) , GTWALL ( 4 .1 ) , GCFNI X ( 20, i ) , 

2 OENTH ( 20 , 1 ) , GW 1 ( 20 , 1 ) , GV J. ( 20 , 1 > , SC ( 3 ) 

D I MENS I ON GAN ( 20 . t > , G AH T ( 20 , 1 ) , OAH 1 ( 20 - 1 ) , GAE X ( 20 , .1. ) , OCPNW ( 4 

1 GPMl-J (41), GPAX ( 4 1 ) , GTNW (41), GT A X (41), GZW ( 4 1 ) , GA*-!GWL ( 4 I. ) , GYV ( 

2 , GWCAR ( 20 , 1 ) , GHCARL < 20 , 1 > , GVCAR i 20 , 1 ) , GARNWL (41), GHTCOE ( 4 1 > 

3 GWMOL ( 20 ) , GMACH ( 20 , 1 ) 

C GEOMETRY DATA OF CELL VERTICES 

DIMENSION G X C ( 2 , 21 ) , GYC( 2,21 ) , GZC(2, 21 ) , GZCELL(42) , GYWALL ( 42 
1 GDYMY ( 42 ) 


1 > : 
20 ) 


) , 
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C U s E R E ■:! L ! I V A L E N C E 3 

EQUIVALENCE ( SC ( 1 ) , RE ().)),( FtfFE, RE ( A ) ) 

DATA MLS Pi N ISP, NRSF*/ 1 , j , .1 / 

DA FA 1. V A R , WAR - CM, VM, il~.HU/ ]. UU^ . u / 

DATA R0A3/8305, 3/ 

1 1 AT A I SURF 1 , I SMFLA/ 1 , 3 / 

C USER PLACES HIS DAI A STATEMENTS HERE. 

BA I A ARcF , URCF : AKcR , L-RCR / 3 • 7 iooE 7 , SUO . , 2 . S283L. 2 , '334 4 ■’ .) . / 

C X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X. X X X X I J S E R SE 0 T 1 1 :• I -J .1. E N I J 3 . 

CX X X X X X X X X X X X X X X X X X X X X XXX XX X X X X XXX X X X X X X STANDARD SECT I Of.,' 2 STAR 1' 
C PLEASE DO HOT ALTER , OR RE -SET, ANY OF THE REMAINING 

0 STATEMENTS OF THIS SECTION. 

IF(SFDATA) 

&CALL RD3FC ( I RN , I NT OR ( 12), LSFDA , MLSP , I SFDA , M I SF , RSFDA , NRSP ) 
ZWTENP=ZWDIST(S1 ) 

ZWDIST<31 ) =0 . 

CALL GRDUTY ( I FT!, I CHAP, I ZED, I MDV-'AP > 

ZWDISTOl ) ==ZM TEMP 

CALL BFCGRD ( I RN , I CHAP , I SUP , I Z ED , I NDVAR > 

I F ( I CHAP .• EG ; -5 ) GO TO J 0 
I F ( I CHAP . LE . 0 . OP , T CHAP . GT . 1 6 ) RETURN 
00 T O ( 1 00 , 200 , 300 , 400 , 500 , 300 , 7 00 , 800 ? 900 , .1. 000 , 1 .1 00 , .1. 200 , 
1300, 1400, 1500, 1300) , I CHAP 
RETURN 

CX XX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X XX X STANDARD SECT I ON 2 ENDS . 
CX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X XXX X X X X USER SEC! I ON 2 STARTS J 


C CHAPTER 0: MODIFY SATLIT DATA, AT START OF EACH I RN 

C I F ( . N 0 T . M A h ! !... 3 T ) R E T U R M 

0 . IF ( I RN , EQ . NRIJM ) DAT F I L- . FALSE . 

C READ SATLIT DATA NAMELIST HERE- 

C CALL l-JR I T 40 ( 4 GHENT ER NAMELIST DATA FOR GROUPS 1 TO. 24 ) 

C READ (20, 01 024) 

C CALL WRIT 40 < 40HEMTER NAMELIST DATA FOR GROUPS 23 TO 42 ) 

C R E A D ( 20,02 5 G 4 2 ) 

10 CONTINUE 

OF I =4. *ATAN< 1 . ) 

C IN NO Z SAT X -RANGE IS SET AS -■ , 01*2. RADIANS 

0 P I 1 0 0 •- G FT * 1 0 0 . 

C OPEN BFCXYZ FILE TO READ GRID VERTICES 


NXP 1 =NX+ 1 
NYP 1 =NY+ 1 
NZF'l =MZ + .l. 

NXNY=NX*NY 
LREC=4*3*NXF 1 *NYP 1 

OPEN ( 23 , P I L.E=" BFCXYZ " , ACCESS* "DIRECT” , BLOCKS I ZE=LREC+4 , 

1 COUNTRY- " RECORD " , STATUS* " OLD " , RECL«LR£C ) 

WRITE (3, 3341 ). 

3541 FORMAT ("FILE BFCXYZ ON LU23 IS OFEM" ) 

I SK I P A=0 
I SK I PB=0 

CALL l-JAL DP ( - 1 , O , N Z , l , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . , 1 . , 1. . , 1 . , VALUE , COEF > 
F MU B* ( 1 . +8 . *FMF B ) / 9 . 

NR I TE ( 3 , 3298 ) FMFB , FMUB , H 1 SAT 

3298 FORMAT ( " FMFB* " , 1 FE .1 2 . 4 , " FMUB " , .1 PE 1 2 . 4 , " H .1 S " , 1 PE .1 2 . 4 ) 

RETURN 
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c CHAPTER 1 : CALLED AT THE START OF EACH TIME STEP. 

C SET " ITT " HERE WHEN T LAST SET NEGATIVE IN BLOCK DATA. 

C "AT I ME + DT ' GIVES THE END TIME OF THE CURRENT T f ME STEF 

C NOT ACCESSED IF STEADY, OR PARABOLIC. 


1.00 CONTINUE 
RETURN 


C CHAPTER 2: CALLED AT THE START OF EACH SWEEP. 

200 COMT I HUE 

1 F ( I SK I PA . GE ■ 1 ) RETURN 
ISKIPA-ISKJFA+1 

C PREPARE WALL TEMPERATURES 

NTHR0=6 

CALL GET 1 D ( Z W 1 . G 7. MODE, NZ > 

Z THRO=G Z NODE ( N THRO ) 

C A L. L G E ' T 1 D « Z G , G Z N 0 D E , N Z > 


CALL T 

MAL 

PC ( N 

Z 

? 

3 

TWhLL 

, G 7. N 

0 

D 

E 7 

z 

THR 

o 

) 






WR I tr ( 

/- 

02* A ) 


T 

3 

:MR, ZT 

HRO 













FORMAT' 

( " W 

ALL 

T 

E 

M 

P ERAT 

ORES 


s 

e r 


I M 

c 

H 2 A* 

(' IS 

wr 

EP 

* 1 

T 4 . " 

.1 .1 F E .! 2 

. 4 > 



















- SETUP 

GE 

OMET 

R 

Y 


CONSTANTS 

- 

- 

~~ 

— 

— 

.... 







DO 33 2 

0 I 

Z=1 : 

N 

Z 

P 

1 














READ ( 2 

3 7 R 

EC = I 

Z 



+ 1 ) < 

( G X L 

( 

I 

X - 

I 

Y ) 7 

T 

Y= 1. ^ r-r 

r'P 1. ) 

,* I 

X 

1 7 H 

! v - 1 ) 

1 






( 

( GYC 

( 

I 

X 7 

I 

Y ) , 

I 

Y — .1 7 N* 

VP 1 > 

, I 

x = 

.1 : M 

X p X ) 







( 

( G Z C 

( 

I 


I 

Y ) , 

I 

Y— 1 , IT 

VP 1 > 

7 I 

X = 

> 1 , N 

IX PI ) 

WR I TE \ 

o ? 3 

003 ) 


( 

G 

VC « 1 7 

I Y ) , 

I 

Y 

= 3. 


NY ) 








FORMAT 

( ■» 

GYC " 

/ 

( 

1 

F'lOEl 

2 . 4 ) 


) 











GZCELL 

( IZ 

) = G Z 

c 

( 

1 

7 .1 ) 














GY WALL 

( IZ 

) ~GY 

c 

/ 

1 

, NYF 1 

) 













GDYMY 

( IZ 

) ~ . 5 


( 


GYC ( 1 

, NVP 

1 


-G 

T 

C ( < 

7 

NY ) > 






THROAT 

-GY 

WALL 

( 

M 

IT 

HRO ) 














MR I TE< 

/-, , ■ -j 

3 23 ) 


( 

G 

ZCELL 

( I Z > 


I 

7 

3. 

. NZ 

4 

.1 ) 







3 1 2 3 F 0 R M A T ( " G 7 C E LL " / < 1 P 1 0 E 1 2 .4) ) 

C CALCULATE WALL INCLINATION ANGLE , 

DO 31.40 12 -1 , NZ 

3 1 40 GANGWL ( I Z ) - At TAN ( < GY WALL ( I Z+ !. ) — GYWALL (12))/ 

1 (01 CELL ( IZ+1 > -GZCELL ( IZ ) ) +1 . E-?> 

WR I TE ( 3 , 3 1 44 ) ( GANGWL ( IZ), I Z-- 1 , NZ ) 

3 1 44 FORMAT ( " WALL AUG" / ( .1. P .1 OE 1 2 .4) ) 

RETURN 


C CHAPTER 3: CALLED AT THE START OF EACH SLAB? 

C MOT ACCESSED IF PARABOLIC, BUT -'STRIDE-' IS. 

300 CONTINUE 
I PELS =50 
I R RLE— 1 50 

G R E L F: A == F LO AT < I SWP- 1 PELS ) /FLOAT ( I RELE- 1 PELS ) 
OREL R A = AM I N 1 ( 1 . , AMAX 1(0., GRELRA ) ) 

RLXF- . 1 + . 4*GRELRA 

RLXPZ=RLXP 

RLXPXY=RLXP 

DTF ALS ( W 1 ) “THROAT «■ ( . 00 1 + 1 00 . *ORELRA ) 

D T F A L S < V 1 ) - D T F A L S ( W 1 > 

RETURN 


CHAPTER 4 5 CALLED AT THE START OF EACH RE-CALCULATION OF 
VARIABLES PI, . . .C4 AT CURRENT SLAB. 1TN0* ITERATION NUMBER. 


400 CONTINUE 
RETURN 


E- 14 



c- 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


c 


CHAPTER 5: GROUND CALLED WHEN SOURCE TERM IS COMPUTED. 
INDVAR GIVES DEPENDENT VARIABLE IN QUESTION IE. LU .,..04, 
TO ADD SOURCE TO DEPENDENT VARIABLE Cl (SAY) FOR IX=IXF, I. XI 
AND I V= I YF » IVL INSERT STATEMENT: 

IF ( INDVAR. EQ. Cl ) 


&CALL ADD ( INDVAR, I XF, IXL., I YF , I YL, TYPE, CM, Vt'i, CVAR, WAR, NY, NX ) 


NOTES ON •ADD-": 

♦SOURCE- ( CVAR ( I Y , I X ) +AMAX 1. ( 0 . 0 , MASFLO > ) * ( WAR ( IV, I X > -PH I ) , 
WHERE •PHI" IS I N— CELL VALUE OF VARIABLE IM QUESTION. 

* - MASFLO ' = CM (IV, IX)* < VM <IY, IX) -F ) , 

WHERE "P" IS THE IN-CELL PRESSURE. 

*FOR INDVAR= Ml, OR =M2, SOURCE ADDED IS "MASFLO"’ ONLY, 

EXCEPT FOR OMEF'HS— . F . $< MASFLO < 0.0 (IE. OUITLOW) WHEN 
CM ( I Y , IX) IS MULTIPLIED BY Rl*Ol (FOR MI) ?< R2»02 (FOR M2 > . 

*BOTH "CVAR" ?•< "CM" ARE MUTLIPL1ED BY CELL-GEOMETRY QUANTITY- 
DICTATED BY SETTING OF "TYPE" ( -CELL , EAST AREA, , - VOLUME) . 

^TYRE-SPECIFIED AREAS ARE CALCULATED AS IF BLOCKAGE ABSENT, 

BUT "VOLUME" WITH ACCOUNT FOR ITS PRESENCE , 

*FOR ALL SOLVED VARIABLES, INCLUDING Ml ( & M2 WHEN ONEPHS=F ) , 


IF " CM > 0.0 CALL "ADD"; FOR Ml g< M2 ALTHOUGH "CVAR" "WAR" 
HAVE NO SIGNIFICANCE THEY MUST BE ENTERED AS ARGUMENTS. 

# " CVAR " , "WAR"', "CM" ?< "VM " MUST BE DIMENSIONED NY, NX. 


500 CONTINUE 

C WRITE (&, 6789) I S UP, I ZED, INDVAR 

067 89 FORMAT ( " IN CH 5 ISW IZ IV ",3I3> 

DO 502 IY-1 , NY 
CVAR ( TV, 1 ) -0, 

502 WAR ( I Y , 1 ) =0 • 

I F ( I N D V A R . N E . C 1 ) G 0 T 0 510 
CALL GET ( C 1. , CPU , MY , NX ) 

CALL GE T ( H2 , G T EMP , NY , NX ) 

CALL 0ET(D1 , GRH, NY , NX > 

JX=l 

DO 508 IY-1, NY 

C CHECK FUEL AND OXIDANT MASS FRACTION 

FUEL = AM AX 1 ( FMFB , CFU ( I Y , I X ) ) 

FMOX-AMAX 1 ( O. ,8.#(CFU<IY, IX)-FMFB) ) 

WAR ( I Y , IX >=FMFB 

‘ CVAR < I Y, I X ) =ARCF*EXP ( -CRCF/GTEMP ( I Y, I X ) ) *GRH ( IV, IX) *#2 
1 *CFU ( I Y , I X > *FMO X / ( FUEL-FMFB+ 1 . E - 1 0 ) 

C SOURCE TERMS FOR REVERSE REACTION 

0 FMPR- 1 , -FMOX-CFU (I Y , I X ) 

C WAR1 ( IY, IX )=FMUB 

C CVAR 1 ( IY, I X ) ARCRSEXP ( -CRCR/GTEMP ( IY, IX) ) #GRH ( IY, IX)**2 

C 1 *FMPR*FMPR / ( F MUD-CFU ( I Y , I X ) + 1 . E- .1 0 > 


SOS 


tj 


CONT I HUE 


CALL 

CALL 


ADD ( L .1 , .1 , NX , 
ADD (Cl , 1 , MX , 


1 - NY, VOLUME, CM, VM, CVAR, WAR, MY, NX ) 

1 - NY , VOLUME , CM , VM , CVAR 1 , WAR 1 , NY , NX > 


RETURN 
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C*#**#**##***#'**-**-*- W A L L F IJ N C T I 0 M S *•»*#•*#«•«• a-*#****!' 
5 .1. 0 I. F ( I 7 . E D . E Q . f i Z > GO 1 0 5 1 5 

I F ( 1 1 IDVAR . HE , H 1 ) GO TO 5 1 7 

C SAVE NEAR WALL DATA 

5 1 5 CALL GET ( W 1 , GW 1 , NY NX ) 

TF< I ZED. L.T. NZ ) OWFLM -~-GWl (NY, 1 ) 

CALL GET « P 1 , GP , MY , MX ) 

GWFP1 ~GF < NY , 1 ) 

CAI. .L GE T ( F 1 H , GP , NY , H X ) 

GWFP 1 M=OP ( MY , .1 ) 

I F ( I ZED . EQ . N7. ) GWFP 1 H“GUFP 1 
CAL L GET ( D L ORH , NY , MX ) 

GWFD1.=GRH(NY, 1 ) 

CALL GET ( D .1 H , GRH , MY , NX ) 

GWFD 1 H-GRH ( NY , 1 ) 

I F < I ZED . EC! . NZ ) GWFD 1 H--GWFD 1 
GPRL-S I GNA ( 24 ) 


C WRITE (6, 6392) GWFP I , OWFFiH, OMFD1 , GWFDtH, GWFW1 

C6392 FORMAT < " P PH D DH Ml IN WALL ” » 1 P5E 12.4) 

5 1. 7 IFH NDVAR . ME . M 1 ) GO TO 520 

C- - WALL. FUNCTIONS FOR Ml 

GDZ - . 5-S ( GZCELL • 1 ZEU+2 ) -GZCELL ( I ZED > ) 

CALL UAL DP ( I ZED, I SWP, NZ, .1 , GDYNY < I ZED) , EMUL AM , GWFD 1 , OWFD.i H . G 
1. 1. . , 1 . , GWFP 1. , GWFP 1 H , GDZ , VALUE , CGEF ) 

CVARv'MY, .1 ) =COEF 
WAR (MY, 1 )=0. 

CALL. ADD ( W 1 , 1 , NX , NY , NY , NORTH , CM, VM , CVAR , WAR , NY , NX ) 

RETURN 

520 T. F ( I NDVAR . ME . H 1 ) GO 7 0 530 

C WALL FUNCTIONS FOR M.1 

CALL GET ( H 1 , GENTH , t I V , MX > 

OD Z «GZCELL < 7 ZED+ 1 ) -GZCELL < I ZED ) 

ODV - . 5, ;>• ( GDYNY ( I ZED • -i-GDYMY ( I ZEU+i ) 1 

CALL WALDP ( I ZED , I SWP , NZ , 2 , GDY , EMULAM , GWFD i , GWFD 1 M , GWFW 1 , 

1 GPRL , S I OH A ( H 1 ) , GWFP 1 , GWFP 1 H, GDZ , VALUE , CGEF ) 

CVAR (NY, 1 ) =COEF 

WAR ( NY , 1 ) “GENTH ( NY , 1 > +GCPNW ( I ZED > * ( OTWALL ( I ZED > -GTNW ( I ZED ) 
CALL. ADD ( H 1,1, NX , NY , NY , NORTH , CM , VM , CVAR , WAR , NY , MX > 

C SAVE HEAT TRANSFER COEFF . FOR PRINTOUT 

I F •: I S W F . E 9 . . L S W E E P ) G H T C 0 E ( I Z E D ) = CO E F 
RETURN 

530 IF( INUVAR.NE.KE) GO TO 540 
C WALL FUNCTIONS FOR KE 


0 0 Z = G 2 C E LL( IZED+-1 ) - G Z C E l . L ( I Z E D ) 

GDY= . 5* ( GDYNY ( I ZED ) +GDYNY ( I ZED+.1 ) ) 

CALL WALDP ( I ZED, I SWP, NZ , 3, GDY, EMULAM, GWFD 1 , GWFD 1 H, GWFW t , 
r GPRL , S I GMA < KE ) , GWFP 1 , GWFP 1 H , GD Z , VALUE , CGEF ) 

CVAR ( NY , 1 ) =COEF 
WAR (NY, 1 ) -VALUE 


CALL ADDUCE, 1 , NX, NY, NY, NORTH, CM, VM, CVAR, WAP, NY, NX ) 
RETURN 

540 IF( I NDVAR. ME. EP) GO TO 550 

GDZ “GZCELL ( I ZED+ 1 ) -GZCELL ( I ZED ) 

GO Y“ . ( GDYNY ( I ZED ) +GDYNY < I ZED+1 ) ) 

CALL WALDP ( I ZED, I SWP , NZ , 4, GDY, EMULAM, GWFD 1 , GWFD 1H, GWFW 1 
1 GPRL , S I GMA ( EP ) , GWFP 1 , GWFP 1 H , GDZ , VALUE , CGEF > 

CVAR ( NY , 1 > “CGEF 
WAR ( NY , l ) “VALUE 


CALL ADD < EP , 1 , MX , NY , MY , NORTH , CM , VM , CVAR , WAR , NY , NX ) 
RETURN 

550 CON T T MUE 
RETURN 




) 
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0 ~ 


c 


CHAPTER 6: CAL LEO AT THE END OF EACH 
CYCLE COMMENCED AT CHAPTER 4. IT MO = 


VAR I ABLE RECA.LCUL AT T. OH 

ITERATION NUMBER. 


600 CONTINUE 
RETURN 


e CM AFTER 


CALLED AT END OF EACH SLAB-WISE CALCULATION, 


700 CONTINUE 

C CALCULATION OF AUXILIARY VARIABLES 

IF ( JSWP. ME. L SWEEP) RETURN 
CALL 'GET ( 1. 6- 0V:|. , NY, MX ) 

CAL L GET ( 2 1 , ON 1 , NY , NX ) 

CALL GET ( P 1 , GP , NY , NX ) 

I F ( I ZED . EQ . 1 ) CALL GET ( F .1 , GP I 7. .1. , NY , N X > 

CALL GET < M2 , GTEMF , NY., NX ) 

CALL GET (III , ORH , MY, MX > 

CALL GET 1. D < ZW 1. , 07.N , MZ •• 

CAL L GRED 1 ( 32 , I Z ED , GAN , NY , MX > 

GARNIIL ( T. ZED ) =GAN C NY , 1. ) 

I F ( I ZED . EG . 1. ) CALL GET ' AH I OH , OAH 1 , NY , NX ) 

I F ( I ZED . EG . NTHRO ) CALL GRED 1 ( 34 , I ZED , OAHT , NY , NX ) 
IF( IZED.EQ.MZ-l ) CALL GRED 1. ( 34 , J ZED , GAEX ,NY,MX> 

I F ( I ZE I.I . TIE. . N THRU ) GU I Li 745 
OATT-0 . 

0FLXT-=0, 

• DO 740 I Y= 1. , NY 
G A T T G A T T + G A H T < I Y , 1 ) #GF 1 1. 00 

GFL X T=OFL X T+QRH ( IV, 1 ) *GAHT ( IY, 1 >*0W1 ( IV, 1 )*0FI100 
740 CONT I HUE 

7 4 5 I F ( I 7 E D .NE.MZ - ). > GO T 0 7 4 9 
0AEXT--0. 

GFLXE-'O, 


DO 746 IY=1,NY 

GAE X T^OAE XT+OAEX ( IY, 1 ) *0P 1 1.00 

GFL X E=*GFL X E+GRH ( IY, IX ) *OAEX ( IY, I X ) *GW1 ( IY, I ) »GPI 1.00 
746 CONTINUE 
749 CONTINUE 

IF ( I ZED. LI . NZ > RETURN 
WRITE (6, 6671 ) 


6 6 7 1 F 0 R M A T ( / / " ##»•#•#•«•*• v- +<■ w ■«■■>!■ * * # «• * * 0 U T P U ' T S i J M T-1 A R Y * * » # # «• v- s -s- # ?- - f * /' ’ 1 ) 
UFANEG-O . 

•J PAP 0S : — O . 

GFA I NJ=0. 

DO 772 I Z — 1 , NTHRO 

772 GPAMEO--OPANEG+GPNW ( I Z ) #GARNWL ( I Z ) *8 I N ( GANGWL !IZ!) #GP 1 1 00 
DO 773 I Z --WTHR'O-i" 1 , N Z - 1 

773 GPAPOS-OPAFOS+GPNW < I Z ) *GARNWL(I Z ') *S1N « GANGWL (1Z) ) *GFI 3. 00 
DO 774 IY=1,NY 

774 ORA I NJ-GF A I M .J+GP I 2 1 ( I Y , 1 ) *GAH 1 ( I Y , 1 ) *GP I 3. 00 
THRUST =GFAPOS+OFANEG+GFA I NU 

G I MPLS- THRUST / < GF LX 1*9 . 8 1 ) 

WRITEI6, 6761 ) GATT, GAEXT , GPAPOS , GFANEG , GFAINJ, GFLXT , GFLXE , 

1. -THRUST, G I MPLS 

676 1 FORMAT ( " THROAT AREA . ", 1. PE 1 2 . 4 , / " E X. I T AREA ” , 1 PE 1 2 . 4 , / 

1 "PRES FORCE +- ”, 1 PE 12.4/'' PRES FORCE -= 'V1PE12.4/ 

1 "PRES FORCE ON INJ=" , 1FE12. 4/"M FLUX THROUGH THROAT" , l PEI. 2. 4/ 

2 "M FLUX THROUGH EXIT " , 1FE12. 4/" THRUST " , 3. PE 1. 2 . 4 

3 /, "SPEC IMF'ULE IN SEC 'M PE 12. 4) 
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Vi 


64 ' 


CAL j._ Nr;. 1 J )J ( Z L i L»/f Jljf Jfcl - NZ ) 

WRITE fo, 6773) 

6 7 7 • *« P Nf-Ti A f ( 1 7 /_ L» /.MU t N A Mm " - yv 

1 "PW PAX TW TNV TAX' 1 ) 

DO 777 IZ=JNNZ 

WR I TE • 6 , 6774 ) I Z - . 5* ( GZCELL ( I 2 > +GZCEL .L < I Z + l > ) , GZNODE < 17.)- 
}. GVWAL L ( I Z ) , GANOWL ( I Z ) *GF I / 1 80 P , GPNW < I Z ) , GPhX ( I Z ) , 

2 OTWALL ( I Z ) , GTNW (17)? GT AX (I Z ) 

COW r I HUE 

4 FORMAT ( 1 3 - 1 F 1 OE 1 2 , 4 ) 

WR T TE ( 6 - 6775) 

6773 FORMAT (/"OUTPUT SUMMARY IN BRITISH UNITS ••/" I Z " ? 5X , " ZG : TH " , 8X : 
1 YN" ,7Xt" PW ( PS I ) PAX ( PS I ) TW ( R ) TNY ( R ) TAX ( R ) ' 

GPS I =6894 . 757 
DO 778 IZ=1»WZ 

778 NR I TE ( 6 * 6 777 ) I Z , 0 7 MODE (17)- GYNALL (17)- GPNW ( 1 1 ) /GPS I , 

1 GPAX (17) / GPS I , GTWALL ( I Z ) * .1 . 8 ? GTNW ( I Z ) * 1 . 8 ? GTAX ( I Z ) * K 8 
6- 7 7 7 F 0 R M A T ( 1 3 - t P l O E 1. 2 , 4 ) 

0 F R I N T H E A 1 ‘ T H A M S P E R C 0 E F 

WR I TE ( 6 - 6301 ) 

6301 FORMAT ( /"HEAT TRANSFER TO THE WALL INFORMATION PRIN TOUT "/ 
i '• I Z ZND COEF COWF*AR QDOT" ) 

DO 779 IZD-NZ 

G ARHT=GHT COE (ID *GARNWL ( I Z ) #GP I 1 00 
7 79 WR I TE ( 6 , 6732 > I Z , 0 Z NODE ( I Z ) - GH T COE ( I Z ) > GARHT ? 
t GARHT# ( OTWAL L (ID -OTMU (17.)) 

6 792 FORMAT ( 1 3 , 1 P7E J 2 . 4 > 

RETURN 


CHAPTER 8s CALLED AT THE END OF EACH SWEEP? 
NOT ACCESSED IF PARABOLIC. 


300 CONTINUE 
RETURN 

C CHAPTER 95 CALLED AT THE END OF EACH TIME STEP? 

C NOT ACCESSED IF PARABOLIC. 

900 CONTINUE 
RETURN 
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C CHAPTER 10: SET PHASE 1 DENSITY HERE WHEN IRH01--1 IN D 

C SET CURRENT -Z "SLAB" DENSITY, D.1 , IF MSLAB-.T. , 

C EG . I F < MSLAB > CALL SET ( D 1,1, NX , 1 , NY , GD 1 , NY , NX 5 . 

C SET NEXT LARGER- Z 'SLAB' DENSITY, 01 H, IF HSLAB~.T. & P 

C . EG . I F < MSLAB ) CAL L SET ( D .1 H , 1 , NX , 1 , NY , GD 1 H , NY , NX ) . 

C SET D(LM(D1.))/DP (IE, D1.DFU FOR UNSTEADY FLOW, 

C EG . IF< MSL AB 5 CALL SET < D 1 DP , 1 , NX , 1 , MY , OD 1 DP , NY , NX ) . 

1000 CONTINUE 

C DENSITY IS CALCULATED FROM EQ OF STATE AS RO--P*W/ ( R* I ) 

CAL L GET (VI, G V , MY , N X ) 

I F ( I ZED . LT . NZ ) CALL GET ( W 1 , GW , NY , MX ) 

CAL L GET (PI, OP , NY , NX ) 

' CALL GET ( C 1 , CFU , NY , NX 1 
CAL L GET ( H 1 , OENTH , NY , NX > 

CALL GET ( H2 , OTEMP , NY , MX ) 

CALL GBIT ( D 1 , GRH , NY , NX ) 

H THOL-14. 0685 
I X - 1 

DO 1010 1. Y-1. ,MY 
I F ( SOL VAR ( C 1 ) ) THEN 
FMOX-S . * ( CFU ( I Y , I X ) -FMF8 ) 

SC ( .1 ) = ( 1 . -FMGX-CFU ( I Y , I X ) ) / 1 8 . 

SC ( 2 > -CFIJ ( I Y , I X >*.5 
SC ( 3 ) =Ff!0 X / 32 . 

WTMOL™ ( SC ( 1 ) # 1 8 . +SC ( 2 ) *2 - +SC ( 3 ) *32 , ) / < SC ( 1 ) +SC ( 2 ) +SC ( 3 > 
END IF 

GWMOL ( I Y ) -WTMOL 
I SUB-0 

I F ( IV. EQ. NY) ISUB=1 

GEK I M- . 5* ( GW ( I Y , I X ) »0W ( T. Y , I X ) +GV ( I Y-I SUB , I X ) *GV ( I Y- 1 SUB 
0HSTA7--H 1 SAT -GEK I N 

I F ( SOLVAR ( H 1 ) ') G H 3 T A T - 0 F. N T H ( I Y , IX) -GEK I N 

CAL L TEMPER ( OHS TAT , G TEMP ( I Y , I X ) , GT , GCFDR , RGAS , SC , 3 ) 

GRM ( I Y , I X ) — GP ( I Y , I. X ) *WTMGL / ( ROAS*GT ) 

GTE MR ( I Y , I X ) -*GT 
1 0 1 0 GCPM IX ( I Y , IX) - G C P D R * R G A S 

CAL L SET ( 2 .1 , NX , 1. , NY , GRH , NY , NX > _ 

CALL SET ( M2 , 1 , NX , 1 , NY , OTEMP , NY , NX ) 

C — SAVE TEMPERATURES AND PRESSURES, AND ALCULATE MACH NO AT I... 
GTNW ( J. Z ED ) =GTEMP ( NY , 1 ) 

‘ G TAX ( I ZED ) = OTEMP (.1,1) 

GCPNW ( T. ZED ) =GCFM I X ( NY , 1 ) 

GPNW ( I ZED ) =GP ( NY , 1 ) 

, GPAX ( I ZED ) -OP (1,1) 

I F ( I SWF .LT. L SWEEP - 1 ) RETURN 
CALL GET ( W2 , GWCAR , MY , N X ) 

IF ( I ZED. LT. NZ ) CALL GET ( W2L , GWCARL , NY , NX ) 

CALL GET ( V2 , GVCAR , NY , NX ) 
l SIJB=0 

DO 1030 I Y ~ 1 , MY 
OWH=GWCAR( I Y , 1 1 

I F ( I Z ED . EQ . N7_ 1 GWH=GWCARL < I Y , l ) 

Gl-JL-GWCARL ( IY, 1 > 

T. F ( I ZED , EQ , 1 ) GWL=GWCAR ( I Y , 1 ) 

GWAV-- . 5* ( GWL+GWH > 

GVAV - . 5* ( GVCAR ( IY, 1 ) +OVCAR ( I Y- 1 SUB , 1 ) ) 

GSOUMD=SQR T ( 1.3 *8305 . * GTE MR ( I Y , 1 )/ GWMOL ( I Y ) ) 

OMAGH ( I Y , 1 ) -SQRT ( GWAV**2+GVAV**2 ) /GSOUND 
.1050 I SUB - 1 ' 

CALL SET ( C2 , 1 , MX , 1 , NY , GMACH , NY , MX ) 

RETURN 


ATA. 

hEAB'-L 


, I X ) ■> 


AS T Si .'EE c- 
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qi ipp 

nirr r mp 


T'WAL 

P 

f; 

( 

NZ 

.. TUALL . 

7. 

MOD 

- 7 r HRO ) 


D I ME 

NS I ON 

7 

NALL 

i 

w 

L. 

t ^ 

ZNOD ( M7. 

) 




COMMQN/MV 


Z TQM 

IP 

i 


"7 i 

- 7 TQW < 2 

, 

) , T 

qm( 27 ) 

- 

ZIP' 

0 AND 

T 

i;»W D 

A 

T 

At 







DATA 

ntqu / 

.•:2 

7/ 










DATA 

ZTQWD 

/ 

-2 . 4 

o 

4 

•“* 

, - 

2 . 1.348, 

- 

i , q 

407 , - 1 , 

1. 

1 

94 1. , . i 

9 

41,. 

■*;< 

■::« 

8 

1 , 

* 9704 , .1 


552 

6 : j , 9 4 0 

1. 

6 . 7 

926 ■, 7 . 

/ 

629 , 

o 

* 

7 

•~J ”! 

3 ? 9 . 7 0 3 

7 

, 1. I 

, 6444. . i 

1 

19. 

4073, 2 

1 

. 348 

o 


2 

• . 

5231 / 





DAT A 

TQW / 1. 


60. t 

1 

r~ 

o 

i\) * 

- 1 5 1 0 . , 

J. 

5 o 0 

, , 14 70. 

1 

830 

. - 820 . 


790 . 

, 

7 

6 

0. 

, 1 450. , 

J 

260 

. : 1060. 


.1 795 , ? 7 65 , ? 7 45 ™ - 730 , , 720 . ? '/ i 5 , - 7 1 0 , / 

C- SORT T WALL FROM INPUT ZTQW AND 'VON ARRAYS SET AS 

DO 510 I Pp~- 1 , NTQW 
5 1 0 Z TQW < I PP ) - Z 7QWD < I PP > + 1 THRO 

1/1 FT I 7 E ( 6 . 678 1 ) ( Z TON ( I RF* > * T ppr: i , NT ON ) 

67 S 1 FORMA T < ! ‘ 7 TON M / ( 1 P5E i 2 * 4 ) > 

ZNODE-ZNOP ( I ) 




:53/17, 4666, 


: j.490 . - 8^0 , 
50. i 8S0 . . 


DA 1"A ABOVE: 


DO 5 1 2 IP1 = 1 7 M l ON 
1. F ( Z NODE . LT . Z '!' QN ( I P 1 ) ) 00 TO 5 i 4 
512 CONTINUE 
514 I PL -IP 1 

:i: PLM--MAXO ( IFL-1 7 J. ) 

TWALL ( 1 ) -TQM ( I PL ) + ( ZTQM \ I PL ) -ZMODE) / ( ZTQM { I PL > -ZTON ( 1 PLM > ■+■ 

1 1 „ E- 1 0 > * ( TON ( J PLM ) -TON ( 1 PL > > 

DO 51.6 I Z 1. -2 ^ MZ 

I F ? Z NOD ( I Z 1 ) . OT . 7 TQM ( I PL. ) ) 00 tq 5 1 7 
516 TWALL ( I Z 1 > --TON ( J PL. ) + ( 7 TON ( I PL ) -ZNOD (III)) / ( Z TON ( I PL > -ZTQM ( IP! 
.1. • l.E-1 0 ) * ( TQM ( I PLM ) -TWO ( I PL > ) 


517 I 7-IZi-l . 

IFF-!: PL ‘ * 

520 IZ=IZ+1 

I F ( I Z . G T „ NZ ) GO tq 530 
Z NODE = Z MOD ( I Z > 

I F ( 7 MODE . GE , 7 TQW < NTQW ) ) GO TO 526 
DO 522 I PP— I PF ? NTQW 
I I I P 1. r ~M 1 NO ( I pp+ 1 , NT ON > 

1 F ( ZMODE . GE - ZTQW ( I PP > , AMD . ZMODE , LT . ZTON (1 1 1 P 1 ) ) GO TO 524 
522 CONTINUE 
524 I FF-M I MO < I PP , NT ON ) 

I PFP=M I NO < I FF+ 1 , NTQW ) 

TWhLL ( I Z ) = TQW « I PFP ) + ( ZTQM ( I PF p ) -ZMODE ) / ( ZTQM < I PFP ) -ZTQM ( I PF } + 
1 l.E-1 O ) * f TQW (IFF) -TQW ( I PFP ) > 

GO TO 520 

526 TWALL ( I Z ) -TQW ( NTQW > 

GO TO 520 
530 CONTINUE 

DO 5555 I Z ■=■“ 1 ? N Z 

MR 1 TE ( 6 , 6784 > IZ, ZNOD ( I Z > ? TWALL ( I Z > 

5555 CONTINUE 

6734 FORMAT < ** IZ Z 'M 3 ■* 1 PE 1 2 ■ 4 - M TWALL-'SJ. PE 1 2 , 4 > 

RETURN 

END 



/ 


to 



SUBROUTINE WALDF ( I 7. 
1 D 


31-IP - HZ , HF- '1 1 1 - DY , AMU - FHO . PM OH , W , R Pl_ 
VALUE * COEF ) 


ri l . 



C THIS SUBROUTINE CALCULATES THE WALL FUNCTIONS P 0R FLOWS 
C FI CANT AXIAL PRESSURE GRADIENTS, FOR REFERENCE SEE? T CE 
C A.M.O. SMITH "A FINITE , . . " , ASME J. BAS , ENG., 1 370, P 523 
C L AUNDER-SPALD I NO "MATHEMATICAL MODELS OF TURBULENCE", AP 


i-.i r h i 
BE.!; I 
, AL; 

1.97 


3 T ON I - 
AND 
0 SET- 


D I MEL 

s i ore 

S V 

PL < 7 

) 7 

SU p L. < 

5 ) , ss ( 5 i 

DATA 

M IIP 

, CA 

P , TA 

UD 

K : TKM 

AX, TKM IN - 5 - .4 , .3 

IJh F A 

nw i. « 

MH ! 

i H V ‘ T: 

, HEP / 1 . , 


■ F 1 RS 

T EM 

FRY 

W J T 

H 

i z™-- 1 

TO SETUP CONS TAM 'VS 

I F < I 7 . 

. GE . 

0 ) 

GO ■} 

rj 

1. 0 0 


01 16 - 

1 1 .. 6 

>- O 





RCAP = 

1 .. / U 

A p 





halo: 

N=. . 

*SG* 

RT ( . 

3 ) 

/CAP 


Art-- 1 . 

+ ! . / 

7 . 





AC= ( 1 

, / y „ 

74 ) 

> - c * ( 2 

. 

A A ) 


ec= < i 

, - :! , 

/ 7 

) / ( J 

„ -f- 

1 . / 7 „ 

) 

CO -4 . 

-ALU 

0 ( 4 

. ) 




i:ap2« 

CAP # 

2 





T2L«0 







RE TUP 

M 






MEN 

SLAB 

EM 

TERE 

D 

SAVE 

3 E I. . E C T E D V A E I A B L . E 3 

I F < I Z 

„EGL 

IZL 

) ’■ j 0 

T 

- j .1 1 0 


I F < 1 7 . 

. EQ , 

1 ) 





I F < I Z 

# p i :i r 

.1 ) 

P L~P 




WAV- „ 

5 * ( W 

+-WL 

) 




DPD2W 

- ( PH 

-P 

) /DZ 




OF- 070 

r ~ ( PH 

-PL 

) / ( 2 

. * 

DZ > 


JF< IZ 

. EH . 

HZ ) 

DPD 

ZG 

= ( p-P 

L ) / DZ 


PL-P 
WL =W 
IZL=IZ 

1 1 0 CONTINUES 

C STAGGERING 

F ( MPH I . ME . MW l ) GO TO 1 20 
PRT~ 1 . 

PRl.~ 1 
PFUN» 1 . 

WP=W 

RHOF= . 5* ( RHO+RHOH ) 

DPD.ZF'-DF'DZW 
GO to 1.21 
120 WP=WAV 
RHOP=RHO 
DF’DZP-DPDZG 

P F U N - 9 . * ( F: R L / P R T - 1 , ) * < F R ' I" / P R I . ) * * . 2 5 
P FUN = P R L. ■*> . 6 6 6 6 6 6 7 

121. PPC--ABS < DP 0 7 . P ) «AMU / < < MP*RHOP > *TIF ) 

RE-RHOF DY^WP / AMU 
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\. 




r 


r 


I F ( MPH I . ME , MW 1 , AMD . MPH I - ME , MH 1 ) GO TO 300 

WALL FUNCTIONS FOR W1 AMD Hi 

PFL-PFC*!7!1 16 

SO “SOR T ( 1 . 4-ppl* 1 1.6) 


AR 

o 

~ 1 l , 6 v- C 

AP- 

V . >r 

( SO - 1 . 

) 


. IF 

/ 

APS ( ARG 

> . L 

T .. 5 

0, ) GO 

T 

0 125 



WRITE ( 

•'j* ■> 

0 l > 

T. Z , HP 

‘HI 

, ARG 

60 1 


FORMAT 

( " j. 

7 . M 

PHI ", 

21 

3 , » BAD ARG IN CALC 



A R Li A M 

I N 1 

( 50 

. ? AM AX 

:L ( 

-50 . , ARG ) ) 

125 E= 


2. +PPL* 

1 1 . 

6+2 

. *SQ) 


EXP (ARG) *.25/11.6 


E»AMhX1 ( 8 . ; E) 

ITERATIVE U+ V + CALCULATION (U+ GUESS=25 ??????? 

VALUE-0. 

CF—AMU / <PFUN*DY) 

S=AC*AMAX 3 ( RE ? 1 . ) ** ( BC— 1 . ) 

}.' F ( RE . LT , 1 32 . 25 ) GO TO 2 1 2 
DO 210 I TR- l » N I TR 
SO- 1 . +30 F: ( j. , +PPC*RE/S > 

SH—SQRT < S > 

L A F - ' A / ( A L. I J- *.? ( 1 , + E * F-! E * 5 H ) — L- L - + A. * * ( -z< ■_• ! — A L LI U ( 5 L! ) ) ) >■ * A 
■:=;:= A MAX 1 ( 1 . E-.1 0» S ) 

210 CONTINUE 


1 PE 1 




CF=CF *RE*S 


■12 

COEF-CF* < 1 

, + P P i„ 

*RE/ 

8 






RETURN 








— 

- WALL FUMC 

TIOMS 

FOR 

‘FT 

<E 




00 

IF (MR HI . ME 

. MKE ) 

GO 

TO 

40 0 





F'F'L-PPC !-vO ). 

16 








SO = 3 ORT( i 

. +PPL 

* 1 .1 . 

6 ) 






ARG— 1. 1 . 6#C 

AP-2 . 

-p- ( C;i7i 

-i 

) 





I F ( APS ( ARG 

) . L T . 

50. ) 

G>! 

j TO 3 

05 




WR I TE ! 

6 1 602 

) T Z 

, Mr HI , AR 

G 


- 

02 

FORMAT 

( " I z 

MPH I 

" 

2 1 3 , " 

BAD 

ARG IN C 

ALC E WALL", 1 PEI 


ARG- AM 

I N 1 ( 5 

0 . ? Al 

MA> 

1 ( -50 

. , AR'G 

) ) 


05 

E- ( 2 . +PPL* 

1 1. . 6+ 

2 . *S 

Q) 

* EXP 

( ARG ) 

*.25/11 

. 6 


E-AMAX 1 (S. 

, p ) 



- 




— 

- ITERATIVE 

LJ+ 

Y + 

CAL 

XULAT 

I ON 

( IJ+ CUES 

•“ “ A 5 •*" •*' ‘ ‘ ) 


S= AC* AM AX 1 

(RE, .1 

. > ** 

( pi: 

: - 1 . ) 





DO 310 ITR 

= 1 , M I 

TR 







SQ-.1 . +SQR T 

( 1 . +F 

FC*Rt/:: 

:> 





SH-SQRT ( S ) 









S--CAP2 / ( AlJ 

GG ( 1 . 

+E*F: 

E#S 

;h ) -cc 

+2. *< 

SQ-ALOG ( 

SQ) ) >**2 


S- AM AX 1 ( 1 .E-10,S> 

310 CONTINUE 

TAU~ S » RHO P * WP * * 2 

TKF I X -AM A X 1 ( TKM I N , TAU/ ( RHOP*TAUOK ) ) 
TKF I X -An I N 1 ( T KM AX - 1 KF I X ) 

VALUE- TKF I X 
COEF “ 1 . E J 0 
RETURN 

WALL FUNCTION FOR ED 

400 VALUE-WAL CON*TKF I X *SQRT ( TKF I X ) / DY 
C0EF=1 , E10 
RETURN 
END 


4 


4 i 
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